]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
A boundary mass matrix does not make much sense in 1D.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 5 Nov 2009 20:57:09 +0000 (20:57 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 5 Nov 2009 20:57:09 +0000 (20:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@20042 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/vectors.templates.h
deal.II/deal.II/source/numerics/matrices.cc

index 266c76ee1f46122046a02ab70aa38b2fce04331e..36745787170c8c825e555dc02d21ad38ae0094f3 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 
 template <class VECTOR, class DH>
-void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>    &mapping,
-                              const DH              &dof,
-                              const Function<DH::space_dimension>   &function,
-                              VECTOR                &vec)
+void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension> &mapping,
+                              const DH                            &dof,
+                              const Function<DH::space_dimension> &function,
+                              VECTOR                              &vec)
 {
   const unsigned int dim=DH::dimension;
 
   Assert (dof.get_fe().n_components() == function.n_components,
          ExcDimensionMismatch(dof.get_fe().n_components(),
                               function.n_components));
-  
+
   const hp::FECollection<DH::dimension,DH::space_dimension> fe (dof.get_fe());
   const unsigned int          n_components = fe.n_components();
   const bool                  fe_is_system = (n_components != 1);
-  
+
   typename DH::active_cell_iterator cell = dof.begin_active(),
                                    endc = dof.end();
 
@@ -88,11 +88,11 @@ void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>
       Assert (unit_support_points[fe_index].size() != 0,
              ExcNonInterpolatingFE());
     }
-  
 
-                                  // Find the support points 
+
+                                  // Find the support points
                                   // on a cell that
-                                  // are multiply mentioned in 
+                                  // are multiply mentioned in
                                   // unit_support_points.
                                   // Mark the first representative
                                   // of each multiply mentioned
@@ -104,7 +104,7 @@ void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>
 
                                   // the following vector collects all dofs i,
                                   // 0<=i<fe.dofs_per_cell, for that
-                                  // unit_support_points[i] 
+                                  // unit_support_points[i]
                                   // is a representative one. i.e.
                                   // the following vector collects all rep dofs.
                                   // the position of a rep dof within this vector
@@ -115,7 +115,7 @@ void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>
   std::vector<std::vector<unsigned int> > dof_to_rep_index_table(fe.size());
 
   std::vector<unsigned int> n_rep_points (fe.size(), 0);
-  
+
   for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
     {
       for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
@@ -128,14 +128,14 @@ void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>
                                           // support points that are placed
                                           // one after the other.
          for (unsigned int j=dofs_of_rep_points[fe_index].size(); j>0; --j)
-           if (unit_support_points[fe_index][i] 
+           if (unit_support_points[fe_index][i]
                == unit_support_points[fe_index][dofs_of_rep_points[fe_index][j-1]])
              {
                dof_to_rep_index_table[fe_index].push_back(j-1);
                representative=false;
                break;
              }
-      
+
          if (representative)
            {
                                               // rep_index=dofs_of_rep_points.size()
@@ -151,12 +151,12 @@ void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>
       Assert(dof_to_rep_index_table[fe_index].size()==fe[fe_index].dofs_per_cell,
             ExcInternalError());
     }
-  
+
   const unsigned int max_rep_points = *std::max_element (n_rep_points.begin(),
                                                         n_rep_points.end());
   std::vector<unsigned int> dofs_on_cell (fe.max_dofs_per_cell());
   std::vector<Point<DH::space_dimension> >  rep_points (max_rep_points);
-  
+
                                   // get space for the values of the
                                   // function at the rep support points.
                                   //
@@ -175,21 +175,21 @@ void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>
                                   // Transformed support points are computed by
                                   // FEValues
   hp::MappingCollection<dim,DH::space_dimension> mapping_collection (mapping);
-  
+
   hp::FEValues<dim, DH::space_dimension> fe_values (mapping_collection,
                               fe, support_quadrature, update_quadrature_points);
-  
+
   for (; cell!=endc; ++cell)
     {
       const unsigned int fe_index = cell->active_fe_index();
-      
+
                                       // for each cell:
                                       // get location of finite element
                                       // support_points
       fe_values.reinit(cell);
       const std::vector<Point<DH::space_dimension> >& support_points =
        fe_values.get_present_fe_values().get_quadrature_points();
-      
+
                                       // pick out the representative
                                       // support points
       rep_points.resize (dofs_of_rep_points[fe_index].size());
@@ -197,7 +197,7 @@ void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>
        rep_points[j] = support_points[dofs_of_rep_points[fe_index][j]];
 
                                       // get indices of the dofs on this cell
-      dofs_on_cell.resize (fe[fe_index].dofs_per_cell);      
+      dofs_on_cell.resize (fe[fe_index].dofs_per_cell);
       cell->get_dof_indices (dofs_on_cell);
 
 
@@ -236,7 +236,7 @@ void VectorTools::interpolate (const Mapping<DH::dimension,DH::space_dimension>
                                           // values to the global
                                           // vector
          for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
-           vec(dofs_on_cell[i]) 
+           vec(dofs_on_cell[i])
              = function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]];
        }
     }
@@ -249,7 +249,7 @@ void VectorTools::interpolate (const DH              &dof,
                               VECTOR                &vec)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  interpolate(StaticMappingQ1<DH::dimension, DH::space_dimension>::mapping, 
+  interpolate(StaticMappingQ1<DH::dimension, DH::space_dimension>::mapping,
              dof, function, vec);
 }
 
@@ -269,27 +269,27 @@ VectorTools::interpolate (const DoFHandler<dim,spacedim>           &dof_1,
 
   std::vector<short unsigned int> touch_count (dof_2.n_dofs(), 0);
   std::vector<unsigned int>       local_dof_indices (dof_2.get_fe().dofs_per_cell);
-  
+
   typename DoFHandler<dim,spacedim>::active_cell_iterator h = dof_1.begin_active();
   typename DoFHandler<dim,spacedim>::active_cell_iterator l = dof_2.begin_active();
   const typename DoFHandler<dim,spacedim>::cell_iterator endh = dof_1.end();
-  
+
   for(; h != endh; ++h, ++l)
   {
     h->get_dof_values(data_1, cell_data_1);
     transfer.vmult(cell_data_2, cell_data_1);
 
     l->get_dof_indices (local_dof_indices);
-  
+
                                   // distribute cell vector
-    for (unsigned int j=0; j<dof_2.get_fe().dofs_per_cell; ++j) 
+    for (unsigned int j=0; j<dof_2.get_fe().dofs_per_cell; ++j)
       {
        data_2(local_dof_indices[j]) += cell_data_2(j);
 
                                         // count, how often we have
                                         // added to this dof
        Assert (touch_count[local_dof_indices[j]] < 255,
-               ExcInternalError());    
+               ExcInternalError());
        ++touch_count[local_dof_indices[j]];
       };
   };
@@ -301,7 +301,7 @@ VectorTools::interpolate (const DoFHandler<dim,spacedim>           &dof_1,
     {
       Assert (touch_count[i] != 0,
              ExcInternalError());
-      
+
       data_2(i) /= touch_count[i];
     };
 }
@@ -315,7 +315,7 @@ namespace internal
 
     template <int dim>
     void
-    interpolate_zero_boundary_values (const dealii::DoFHandler<dim>         &dof_handler,
+    interpolate_zero_boundary_values (const dealii::DoFHandler<dim> &dof_handler,
                                       std::map<unsigned int,double> &boundary_values)
     {
                                        // we only need to find the
@@ -337,8 +337,8 @@ namespace internal
 
 //codimension 1
     void
-      interpolate_zero_boundary_values (const dealii::DoFHandler<1,2>         &dof_handler,
-                                      std::map<unsigned int,double> &boundary_values)
+      interpolate_zero_boundary_values (const dealii::DoFHandler<1,2> &dof_handler,
+                                       std::map<unsigned int,double> &boundary_values)
     {
                                        // we only need to find the
                                        // left-most and right-most
@@ -358,11 +358,11 @@ namespace internal
 
 
 #else
-    
+
     template <int dim, int spacedim>
     void
-    interpolate_zero_boundary_values (const dealii::DoFHandler<dim,spacedim>       &dof_handler,
-                                      std::map<unsigned int,double> &boundary_values)
+    interpolate_zero_boundary_values (const dealii::DoFHandler<dim,spacedim> &dof_handler,
+                                      std::map<unsigned int,double>          &boundary_values)
     {
       const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
 
@@ -405,7 +405,7 @@ namespace internal
              boundary_values[face_dof_indices[i]] = 0.;
          }
     }
-    
+
 #endif
   }
 }
@@ -413,15 +413,15 @@ namespace internal
 
 
 template <int dim, class VECTOR, int spacedim>
-void VectorTools::project (const Mapping<dim, spacedim>       &mapping,
-                          const DoFHandler<dim,spacedim>    &dof,
-                          const ConstraintMatrix   &constraints,
-                          const Quadrature<dim>    &quadrature,
-                          const Function<spacedim>      &function,
-                          VECTOR                   &vec_result,
-                          const bool                enforce_zero_boundary,
-                          const Quadrature<dim-1>  &q_boundary,
-                          const bool                project_to_boundary_first)
+void VectorTools::project (const Mapping<dim,spacedim>     &mapping,
+                          const DoFHandler<dim,spacedim>  &dof,
+                          const ConstraintMatrix          &constraints,
+                          const Quadrature<dim>           &quadrature,
+                          const Function<spacedim>        &function,
+                          VECTOR                          &vec_result,
+                          const bool                       enforce_zero_boundary,
+                          const Quadrature<dim-1>         &q_boundary,
+                          const bool                       project_to_boundary_first)
 {
   Assert (dof.get_fe().n_components() == function.n_components,
          ExcDimensionMismatch(dof.get_fe().n_components(),
@@ -429,18 +429,18 @@ void VectorTools::project (const Mapping<dim, spacedim>       &mapping,
 
   Assert (vec_result.size() == dof.n_dofs(),
           ExcDimensionMismatch (vec_result.size(), dof.n_dofs()));
-  
+
                                   // make up boundary values
   std::map<unsigned int,double> boundary_values;
-  
-  if (enforce_zero_boundary == true) 
+
+  if (enforce_zero_boundary == true)
                                     // no need to project boundary
                                     // values, but enforce
                                     // homogeneous boundary values
                                     // anyway
     internal::VectorTools::
       interpolate_zero_boundary_values (dof, boundary_values);
-  
+
   else
                                     // no homogeneous boundary values
     if (project_to_boundary_first == true)
@@ -500,7 +500,7 @@ void VectorTools::project (const Mapping<dim, spacedim>       &mapping,
   prec.initialize(mass_matrix, 1.2);
                                   // solve
   cg.solve (mass_matrix, vec, tmp, prec);
-  
+
                                   // distribute solution
   constraints.distribute (vec);
 
@@ -517,14 +517,14 @@ void VectorTools::project (const Mapping<dim, spacedim>       &mapping,
 
 
 template <int dim, class VECTOR, int spacedim>
-void VectorTools::project (const DoFHandler<dim,spacedim>    &dof,
-                          const ConstraintMatrix   &constraints,
-                          const Quadrature<dim>    &quadrature,
-                          const Function<spacedim>      &function,
-                          VECTOR                   &vec,
-                          const bool                enforce_zero_boundary,
-                          const Quadrature<dim-1>  &q_boundary,
-                          const bool                project_to_boundary_first)
+void VectorTools::project (const DoFHandler<dim,spacedim> &dof,
+                          const ConstraintMatrix         &constraints,
+                          const Quadrature<dim>          &quadrature,
+                          const Function<spacedim>       &function,
+                          VECTOR                         &vec,
+                          const bool                      enforce_zero_boundary,
+                          const Quadrature<dim-1>        &q_boundary,
+                          const bool                      project_to_boundary_first)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   project(StaticMappingQ1<dim,spacedim>::mapping, dof, constraints, quadrature, function, vec,
@@ -535,11 +535,12 @@ void VectorTools::project (const DoFHandler<dim,spacedim>    &dof,
 
 
 template <int dim, int spacedim>
-void VectorTools::create_right_hand_side (const Mapping<dim, spacedim>    &mapping,
-                                         const DoFHandler<dim,spacedim> &dof_handler,
-                                         const Quadrature<dim> &quadrature,
-                                         const Function<spacedim>   &rhs_function,
-                                         Vector<double>        &rhs_vector)
+void VectorTools::create_right_hand_side
+  (const Mapping<dim,spacedim>    &mapping,
+   const DoFHandler<dim,spacedim> &dof_handler,
+   const Quadrature<dim>          &quadrature,
+   const Function<spacedim>       &rhs_function,
+   Vector<double>                 &rhs_vector)
 {
   const FiniteElement<dim,spacedim> &fe  = dof_handler.get_fe();
   Assert (fe.n_components() == rhs_function.n_components,
@@ -547,7 +548,7 @@ void VectorTools::create_right_hand_side (const Mapping<dim, spacedim>    &mappi
   Assert (rhs_vector.size() == dof_handler.n_dofs(),
          ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
   rhs_vector = 0;
-  
+
   UpdateFlags update_flags = UpdateFlags(update_values   |
                                         update_quadrature_points |
                                         update_JxW_values);
@@ -556,7 +557,7 @@ void VectorTools::create_right_hand_side (const Mapping<dim, spacedim>    &mappi
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                     n_q_points    = fe_values.n_quadrature_points,
                     n_components  = fe.n_components();
-  
+
   std::vector<unsigned int> dofs (dofs_per_cell);
   Vector<double> cell_vector (dofs_per_cell);
 
@@ -567,42 +568,42 @@ void VectorTools::create_right_hand_side (const Mapping<dim, spacedim>    &mappi
   if (n_components==1)
     {
       std::vector<double> rhs_values(n_q_points);
-      
-      for (; cell!=endc; ++cell) 
+
+      for (; cell!=endc; ++cell)
        {
          fe_values.reinit(cell);
-         
+
          const std::vector<double> &weights   = fe_values.get_JxW_values ();
          rhs_function.value_list (fe_values.get_quadrature_points(),
                                   rhs_values);
-         
+
          cell_vector = 0;
          for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i) 
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
              cell_vector(i) += rhs_values[point] *
                                fe_values.shape_value(i,point) *
                                weights[point];
-       
+
          cell->get_dof_indices (dofs);
-         
+
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            rhs_vector(dofs[i]) += cell_vector(i);
        }
-      
+
     }
   else
     {
       std::vector<Vector<double> > rhs_values(n_q_points,
                                              Vector<double>(n_components));
-      
-      for (; cell!=endc; ++cell) 
+
+      for (; cell!=endc; ++cell)
        {
          fe_values.reinit(cell);
-         
+
          const std::vector<double> &weights   = fe_values.get_JxW_values ();
          rhs_function.vector_value_list (fe_values.get_quadrature_points(),
                                          rhs_values);
-             
+
          cell_vector = 0;
                                           // Use the faster code if the
                                           // FiniteElement is primitive
@@ -613,7 +614,7 @@ void VectorTools::create_right_hand_side (const Mapping<dim, spacedim>    &mappi
                  {
                    const unsigned int component
                      = fe.system_to_component_index(i).first;
-                   
+
                    cell_vector(i) += rhs_values[point](component) *
                                      fe_values.shape_value(i,point) *
                                      weights[point];
@@ -634,9 +635,9 @@ void VectorTools::create_right_hand_side (const Mapping<dim, spacedim>    &mappi
                                          weights[point];
                      }
            }
-      
+
          cell->get_dof_indices (dofs);
-         
+
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            rhs_vector(dofs[i]) += cell_vector(i);
        }
@@ -646,25 +647,27 @@ void VectorTools::create_right_hand_side (const Mapping<dim, spacedim>    &mappi
 
 
 template <int dim, int spacedim>
-void VectorTools::create_right_hand_side (const DoFHandler<dim,spacedim>    &dof_handler,
-                                         const Quadrature<dim>    &quadrature,
-                                         const Function<spacedim>      &rhs_function,
-                                         Vector<double>           &rhs_vector)
+void VectorTools::create_right_hand_side
+  (const DoFHandler<dim,spacedim> &dof_handler,
+   const Quadrature<dim>          &quadrature,
+   const Function<spacedim>       &rhs_function,
+   Vector<double>                 &rhs_vector)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  create_right_hand_side(StaticMappingQ1<dim,spacedim>::mapping, dof_handler, quadrature,
-                        rhs_function, rhs_vector);
+  create_right_hand_side(StaticMappingQ1<dim,spacedim>::mapping, dof_handler,
+                        quadrature, rhs_function, rhs_vector);
 }
 
 
 
 
 template <int dim, int spacedim>
-void VectorTools::create_right_hand_side (const hp::MappingCollection<dim,spacedim>    &mapping,
-                                         const hp::DoFHandler<dim,spacedim> &dof_handler,
-                                         const hp::QCollection<dim> &quadrature,
-                                         const Function<spacedim>   &rhs_function,
-                                         Vector<double>        &rhs_vector)
+void VectorTools::create_right_hand_side
+  (const hp::MappingCollection<dim,spacedim> &mapping,
+   const hp::DoFHandler<dim,spacedim>        &dof_handler,
+   const hp::QCollection<dim>                &quadrature,
+   const Function<spacedim>                  &rhs_function,
+   Vector<double>                            &rhs_vector)
 {
   const hp::FECollection<dim,spacedim> &fe  = dof_handler.get_fe();
   Assert (fe.n_components() == rhs_function.n_components,
@@ -672,14 +675,14 @@ void VectorTools::create_right_hand_side (const hp::MappingCollection<dim,spaced
   Assert (rhs_vector.size() == dof_handler.n_dofs(),
          ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
   rhs_vector = 0;
-  
+
   UpdateFlags update_flags = UpdateFlags(update_values   |
                                         update_quadrature_points |
                                         update_JxW_values);
   hp::FEValues<dim,spacedim> x_fe_values (mapping, fe, quadrature, update_flags);
 
   const unsigned int n_components  = fe.n_components();
-  
+
   std::vector<unsigned int> dofs (fe.max_dofs_per_cell());
   Vector<double> cell_vector (fe.max_dofs_per_cell());
 
@@ -690,60 +693,62 @@ void VectorTools::create_right_hand_side (const hp::MappingCollection<dim,spaced
   if (n_components==1)
     {
       std::vector<double> rhs_values;
-      
-      for (; cell!=endc; ++cell) 
+
+      for (; cell!=endc; ++cell)
        {
          x_fe_values.reinit(cell);
 
-         const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
-         
+         const FEValues<dim,spacedim> &fe_values =
+           x_fe_values.get_present_fe_values();
+
          const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                             n_q_points    = fe_values.n_quadrature_points;
          rhs_values.resize (n_q_points);
          dofs.resize (dofs_per_cell);
          cell_vector.reinit (dofs_per_cell);
-         
+
          const std::vector<double> &weights   = fe_values.get_JxW_values ();
          rhs_function.value_list (fe_values.get_quadrature_points(),
                                   rhs_values);
-         
+
          cell_vector = 0;
          for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i) 
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
              cell_vector(i) += rhs_values[point] *
                                fe_values.shape_value(i,point) *
                                weights[point];
-       
+
          cell->get_dof_indices (dofs);
-         
+
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            rhs_vector(dofs[i]) += cell_vector(i);
        }
-      
+
     }
   else
     {
       std::vector<Vector<double> > rhs_values;
-      
+
       for (; cell!=endc; ++cell)
        {
          x_fe_values.reinit(cell);
 
-         const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
-         
+         const FEValues<dim,spacedim> &fe_values =
+           x_fe_values.get_present_fe_values();
+
          const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                             n_q_points    = fe_values.n_quadrature_points;
          rhs_values.resize (n_q_points,
                             Vector<double>(n_components));
          dofs.resize (dofs_per_cell);
          cell_vector.reinit (dofs_per_cell);
-             
+
          const std::vector<double> &weights   = fe_values.get_JxW_values ();
          rhs_function.vector_value_list (fe_values.get_quadrature_points(),
                                          rhs_values);
-             
+
          cell_vector = 0;
-         
+
                                           // Use the faster code if the
                                           // FiniteElement is primitive
          if (cell->get_fe().is_primitive ())
@@ -753,7 +758,7 @@ void VectorTools::create_right_hand_side (const hp::MappingCollection<dim,spaced
                  {
                    const unsigned int component
                      = cell->get_fe().system_to_component_index(i).first;
-                   
+
                    cell_vector(i) += rhs_values[point](component) *
                                      fe_values.shape_value(i,point) *
                                      weights[point];
@@ -773,9 +778,9 @@ void VectorTools::create_right_hand_side (const hp::MappingCollection<dim,spaced
                                          weights[point];
                      }
            }
-             
+
          cell->get_dof_indices (dofs);
-         
+
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            rhs_vector(dofs[i]) += cell_vector(i);
        }
@@ -785,13 +790,14 @@ void VectorTools::create_right_hand_side (const hp::MappingCollection<dim,spaced
 
 
 template <int dim, int spacedim>
-void VectorTools::create_right_hand_side (const hp::DoFHandler<dim,spacedim>    &dof_handler,
-                                         const hp::QCollection<dim>    &quadrature,
-                                         const Function<spacedim>      &rhs_function,
-                                         Vector<double>           &rhs_vector)
+void VectorTools::create_right_hand_side
+  (const hp::DoFHandler<dim,spacedim>  &dof_handler,
+   const hp::QCollection<dim>          &quadrature,
+   const Function<spacedim>            &rhs_function,
+   Vector<double>                      &rhs_vector)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  create_right_hand_side(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, 
+  create_right_hand_side(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
                         dof_handler, quadrature,
                         rhs_function, rhs_vector);
 }
@@ -800,16 +806,17 @@ void VectorTools::create_right_hand_side (const hp::DoFHandler<dim,spacedim>
 
 
 template <int dim, int spacedim>
-void VectorTools::create_point_source_vector (const Mapping<dim, spacedim>       &mapping,
-                                              const DoFHandler<dim,spacedim>    &dof_handler,
-                                              const Point<spacedim>         &p,
-                                              Vector<double>           &rhs_vector)
+void VectorTools::create_point_source_vector
+  (const Mapping<dim,spacedim>    &mapping,
+   const DoFHandler<dim,spacedim> &dof_handler,
+   const Point<spacedim>          &p,
+   Vector<double>                 &rhs_vector)
 {
    Assert (rhs_vector.size() == dof_handler.n_dofs(),
            ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
    Assert (dof_handler.get_fe().n_components() == 1,
           ExcMessage ("This function only works for scalar finite elements"));
-   
+
    rhs_vector = 0;
 
    std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
@@ -834,9 +841,10 @@ void VectorTools::create_point_source_vector (const Mapping<dim, spacedim>
 
 
 template <int dim, int spacedim>
-void VectorTools::create_point_source_vector (const DoFHandler<dim,spacedim>    &dof_handler,
-                                              const Point<spacedim>         &p,
-                                              Vector<double>           &rhs_vector)
+void VectorTools::create_point_source_vector
+  (const DoFHandler<dim,spacedim> &dof_handler,
+   const Point<spacedim>          &p,
+   Vector<double>                 &rhs_vector)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   create_point_source_vector(StaticMappingQ1<dim,spacedim>::mapping, dof_handler,
@@ -845,16 +853,17 @@ void VectorTools::create_point_source_vector (const DoFHandler<dim,spacedim>
 
 
 template <int dim, int spacedim>
-void VectorTools::create_point_source_vector (const hp::MappingCollection<dim,spacedim>       &mapping,
-                                              const hp::DoFHandler<dim,spacedim>    &dof_handler,
-                                              const Point<spacedim>         &p,
-                                              Vector<double>           &rhs_vector)
+void VectorTools::create_point_source_vector
+  (const hp::MappingCollection<dim,spacedim> &mapping,
+   const hp::DoFHandler<dim,spacedim>        &dof_handler,
+   const Point<spacedim>                     &p,
+   Vector<double>                            &rhs_vector)
 {
    Assert (rhs_vector.size() == dof_handler.n_dofs(),
            ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
    Assert (dof_handler.get_fe().n_components() == 1,
           ExcMessage ("This function only works for scalar finite elements"));
-   
+
    rhs_vector = 0;
 
    std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
@@ -879,9 +888,10 @@ void VectorTools::create_point_source_vector (const hp::MappingCollection<dim,sp
 
 
 template <int dim, int spacedim>
-void VectorTools::create_point_source_vector (const hp::DoFHandler<dim,spacedim>    &dof_handler,
-                                              const Point<spacedim>         &p,
-                                              Vector<double>           &rhs_vector)
+void VectorTools::create_point_source_vector
+  (const hp::DoFHandler<dim,spacedim> &dof_handler,
+   const Point<spacedim>              &p,
+   Vector<double>                     &rhs_vector)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   create_point_source_vector(hp::StaticMappingQ1<dim>::mapping_collection,
@@ -894,21 +904,22 @@ void VectorTools::create_point_source_vector (const hp::DoFHandler<dim,spacedim>
 
 template <int dim, int spacedim>
 void
-VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>      &mapping,
-                                             const DoFHandler<dim,spacedim>   &dof_handler,
-                                             const Quadrature<dim-1> &quadrature,
-                                             const Function<spacedim>     &rhs_function,
-                                             Vector<double>          &rhs_vector,
-                                             const std::set<unsigned char> &boundary_indicators)
+VectorTools::create_boundary_right_hand_side
+  (const Mapping<dim, spacedim>   &mapping,
+   const DoFHandler<dim,spacedim> &dof_handler,
+   const Quadrature<dim-1>        &quadrature,
+   const Function<spacedim>       &rhs_function,
+   Vector<double>                 &rhs_vector,
+   const std::set<unsigned char>  &boundary_indicators)
 {
   const FiniteElement<dim> &fe  = dof_handler.get_fe();
   Assert (fe.n_components() == rhs_function.n_components,
          ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
   Assert (rhs_vector.size() == dof_handler.n_dofs(),
          ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
-  
+
   rhs_vector = 0;
-  
+
   UpdateFlags update_flags = UpdateFlags(update_values   |
                                         update_quadrature_points |
                                         update_JxW_values);
@@ -917,17 +928,18 @@ VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                     n_q_points    = fe_values.n_quadrature_points,
                     n_components  = fe.n_components();
-  
+
   std::vector<unsigned int> dofs (dofs_per_cell);
   Vector<double> cell_vector (dofs_per_cell);
 
-  typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof_handler.begin_active(),
-                                                endc = dof_handler.end();
+  typename DoFHandler<dim,spacedim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
 
   if (n_components==1)
     {
       std::vector<double> rhs_values(n_q_points);
-      
+
       for (; cell!=endc; ++cell)
        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
          if (cell->face(face)->at_boundary () &&
@@ -936,19 +948,19 @@ VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>
               boundary_indicators.end()))
            {
              fe_values.reinit(cell, face);
-         
+
              const std::vector<double> &weights   = fe_values.get_JxW_values ();
              rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
-             
+
              cell_vector = 0;
              for (unsigned int point=0; point<n_q_points; ++point)
-               for (unsigned int i=0; i<dofs_per_cell; ++i) 
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
                  cell_vector(i) += rhs_values[point] *
                                    fe_values.shape_value(i,point) *
                                    weights[point];
-       
+
              cell->get_dof_indices (dofs);
-         
+
              for (unsigned int i=0; i<dofs_per_cell; ++i)
                rhs_vector(dofs[i]) += cell_vector(i);
            }
@@ -956,8 +968,8 @@ VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>
   else
     {
       std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-      
-      for (; cell!=endc; ++cell) 
+
+      for (; cell!=endc; ++cell)
        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
          if (cell->face(face)->at_boundary () &&
              (boundary_indicators.find (cell->face(face)->boundary_indicator())
@@ -965,22 +977,22 @@ VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>
               boundary_indicators.end()))
            {
              fe_values.reinit(cell, face);
-             
+
              const std::vector<double> &weights   = fe_values.get_JxW_values ();
              rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-             
+
              cell_vector = 0;
-             
+
                                               // Use the faster code if the
                                               // FiniteElement is primitive
              if (fe.is_primitive ())
-               {                 
+               {
                  for (unsigned int point=0; point<n_q_points; ++point)
                    for (unsigned int i=0; i<dofs_per_cell; ++i)
                      {
                        const unsigned int component
                          = fe.system_to_component_index(i).first;
-                       
+
                        cell_vector(i) += rhs_values[point](component) *
                                          fe_values.shape_value(i,point) *
                                          weights[point];
@@ -1002,9 +1014,9 @@ VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>
                              weights[point];
                          }
                }
-                 
+
              cell->get_dof_indices (dofs);
-             
+
              for (unsigned int i=0; i<dofs_per_cell; ++i)
                rhs_vector(dofs[i]) += cell_vector(i);
            }
@@ -1016,12 +1028,12 @@ VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>
 // Implementation for 1D
 template <int dim, int spacedim>
 void
-VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>    &,
+VectorTools::create_boundary_right_hand_side (const Mapping<dim,spacedim>    &,
                                              const DoFHandler<dim,spacedim> &,
-                                             const Quadrature<dim-1> &,
-                                             const Function<spacedim>   &,
-                                             Vector<double>      &,
-                                             const std::set<unsigned char> &)
+                                             const Quadrature<dim-1>        &,
+                                             const Function<spacedim>       &,
+                                             Vector<double>                 &,
+                                             const std::set<unsigned char>  &)
 {
   Assert (false, ExcImpossibleInDim(dim));
 }
@@ -1030,11 +1042,12 @@ VectorTools::create_boundary_right_hand_side (const Mapping<dim, spacedim>    &,
 
 template <int dim, int spacedim>
 void
-VectorTools::create_boundary_right_hand_side (const DoFHandler<dim,spacedim>   &dof_handler,
-                                             const Quadrature<dim-1> &quadrature,
-                                             const Function<spacedim>     &rhs_function,
-                                             Vector<double>          &rhs_vector,
-                                             const std::set<unsigned char> &boundary_indicators)
+VectorTools::create_boundary_right_hand_side
+  (const DoFHandler<dim,spacedim> &dof_handler,
+   const Quadrature<dim-1>        &quadrature,
+   const Function<spacedim>       &rhs_function,
+   Vector<double>                 &rhs_vector,
+   const std::set<unsigned char>  &boundary_indicators)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
 
@@ -1050,28 +1063,29 @@ VectorTools::create_boundary_right_hand_side (const DoFHandler<dim,spacedim>   &
 
 template <int dim, int spacedim>
 void
-VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,spacedim>      &mapping,
-                                             const hp::DoFHandler<dim,spacedim>   &dof_handler,
-                                             const hp::QCollection<dim-1> &quadrature,
-                                             const Function<spacedim>     &rhs_function,
-                                             Vector<double>          &rhs_vector,
-                                             const std::set<unsigned char> &boundary_indicators)
+VectorTools::create_boundary_right_hand_side
+  (const hp::MappingCollection<dim,spacedim> &mapping,
+   const hp::DoFHandler<dim,spacedim>        &dof_handler,
+   const hp::QCollection<dim-1>              &quadrature,
+   const Function<spacedim>                  &rhs_function,
+   Vector<double>                            &rhs_vector,
+   const std::set<unsigned char>             &boundary_indicators)
 {
   const hp::FECollection<dim> &fe  = dof_handler.get_fe();
   Assert (fe.n_components() == rhs_function.n_components,
          ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
   Assert (rhs_vector.size() == dof_handler.n_dofs(),
          ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
-  
+
   rhs_vector = 0;
-  
+
   UpdateFlags update_flags = UpdateFlags(update_values   |
                                         update_quadrature_points |
                                         update_JxW_values);
   hp::FEFaceValues<dim> x_fe_values (mapping, fe, quadrature, update_flags);
 
   const unsigned int n_components  = fe.n_components();
-  
+
   std::vector<unsigned int> dofs (fe.max_dofs_per_cell());
   Vector<double> cell_vector (fe.max_dofs_per_cell());
 
@@ -1082,7 +1096,7 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,sp
   if (n_components==1)
     {
       std::vector<double> rhs_values;
-      
+
       for (; cell!=endc; ++cell)
        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
          if (cell->face(face)->at_boundary () &&
@@ -1097,19 +1111,19 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,sp
              const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                                 n_q_points    = fe_values.n_quadrature_points;
              rhs_values.resize (n_q_points);
-             
+
              const std::vector<double> &weights   = fe_values.get_JxW_values ();
              rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
-             
+
              cell_vector = 0;
              for (unsigned int point=0; point<n_q_points; ++point)
-               for (unsigned int i=0; i<dofs_per_cell; ++i) 
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
                  cell_vector(i) += rhs_values[point] *
                                    fe_values.shape_value(i,point) *
                                    weights[point];
-       
+
              cell->get_dof_indices (dofs);
-         
+
              for (unsigned int i=0; i<dofs_per_cell; ++i)
                rhs_vector(dofs[i]) += cell_vector(i);
            }
@@ -1117,8 +1131,8 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,sp
   else
     {
       std::vector<Vector<double> > rhs_values;
-      
-      for (; cell!=endc; ++cell) 
+
+      for (; cell!=endc; ++cell)
        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
          if (cell->face(face)->at_boundary () &&
              (boundary_indicators.find (cell->face(face)->boundary_indicator())
@@ -1132,22 +1146,22 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,sp
              const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                                 n_q_points    = fe_values.n_quadrature_points;
              rhs_values.resize (n_q_points, Vector<double>(n_components));
-             
+
              const std::vector<double> &weights   = fe_values.get_JxW_values ();
              rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-             
+
              cell_vector = 0;
-             
+
                                               // Use the faster code if the
                                               // FiniteElement is primitive
              if (cell->get_fe().is_primitive ())
-               {                 
+               {
                  for (unsigned int point=0; point<n_q_points; ++point)
                    for (unsigned int i=0; i<dofs_per_cell; ++i)
                      {
                        const unsigned int component
                          = cell->get_fe().system_to_component_index(i).first;
-                       
+
                        cell_vector(i) += rhs_values[point](component) *
                                          fe_values.shape_value(i,point) *
                                          weights[point];
@@ -1169,9 +1183,9 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,sp
                              weights[point];
                          }
                }
-                 
+
              cell->get_dof_indices (dofs);
-             
+
              for (unsigned int i=0; i<dofs_per_cell; ++i)
                rhs_vector(dofs[i]) += cell_vector(i);
            }
@@ -1183,7 +1197,7 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,sp
 // Implementation for 1D
 template <int dim, int spacedim>
 void
-VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,spacedim>    &,
+VectorTools::create_boundary_right_hand_side (const hp::MappingCollection<dim,spacedim> &,
                                              const hp::DoFHandler<dim,spacedim> &,
                                              const hp::QCollection<dim-1> &,
                                              const Function<spacedim>   &,
@@ -1240,7 +1254,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
          ExcMessage ("The number of components in the mask has to be either "
                      "zero or equal to the number of components in the finite "
                      "element."));
-    
+
                                   // check whether boundary values at
                                   // the left or right boundary of
                                   // the line are
@@ -1251,7 +1265,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
                                   // right.
   const unsigned int direction = boundary_component;
   Assert (direction < 2, ExcInvalidBoundaryIndicator());
-  
+
                                   // first find the outermost active
                                   // cell by first traversing the coarse
                                   // grid to its end and then going
@@ -1259,7 +1273,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
   typename DH::cell_iterator outermost_cell = dof.begin(0);
   while (outermost_cell->neighbor(direction).state() == IteratorState::valid)
     outermost_cell = outermost_cell->neighbor(direction);
-  
+
   while (outermost_cell->has_children())
     outermost_cell = outermost_cell->child(direction);
 
@@ -1296,7 +1310,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
   else
     boundary_function.vector_value (outermost_cell->vertex(direction),
                                    function_values);
-  
+
   for (unsigned int i=0; i<fe.dofs_per_vertex; ++i)
     if (component_mask[fe.face_system_to_component_index(i).first])
       boundary_values[outermost_cell->vertex_dof_index(direction,i)]
@@ -1353,7 +1367,7 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
                                   // immediately
   if (function_map.size() == 0)
     return;
-  
+
   Assert (function_map.find(255) == function_map.end(),
          ExcInvalidBoundaryIndicator());
 
@@ -1403,7 +1417,7 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
   for (unsigned int f=0; f<finite_elements.size(); ++f)
     {
       const FiniteElement<dim> &fe = finite_elements[f];
-      
+
                                       // generate a quadrature rule
                                       // on the face from the unit
                                       // support points. this will be
@@ -1460,9 +1474,9 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
              if (component_mask[fe.face_system_to_component_index(i).first]
                  == true)
                unit_support_points[i] = fe.unit_face_support_point(i);
-       
+
          q_collection.push_back (Quadrature<dim-1>(unit_support_points));
-        }    
+        }
     }
                                   // now that we have a q_collection
                                   // object with all the right
@@ -1473,7 +1487,7 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
   hp::MappingCollection<dim> mapping_collection (mapping);
   hp::FEFaceValues<dim> x_fe_values (mapping_collection, finite_elements, q_collection,
                                     update_quadrature_points);
-  
+
   typename DH::active_cell_iterator cell = dof.begin_active(),
                                    endc = dof.end();
   for (; cell!=endc; ++cell)
@@ -1505,10 +1519,10 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
                                    "values that correspond to primitive (scalar) base "
                                    "elements"));
          }
-       
+
        typename DH::face_iterator face = cell->face(face_no);
        const unsigned char boundary_component = face->boundary_indicator();
-       if (function_map.find(boundary_component) != function_map.end()) 
+       if (function_map.find(boundary_component) != function_map.end())
          {
                                             // face is of the right component
            x_fe_values.reinit(cell, face_no);
@@ -1521,7 +1535,7 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
            face->get_dof_indices (face_dofs, cell->active_fe_index());
            const std::vector<Point<DH::space_dimension> > &dof_locations
               = fe_values.get_quadrature_points ();
-           
+
            if (fe_is_system)
              {
                                                  // resize
@@ -1535,10 +1549,10 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
                                             Vector<double>(fe.n_components()));
                 else
                   dof_values_system.resize (fe.dofs_per_face);
-                
+
                function_map.find(boundary_component)->second
                   ->vector_value_list (dof_locations, dof_values_system);
-               
+
                                                 // enter those dofs
                                                 // into the list that
                                                 // match the
@@ -1646,10 +1660,10 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
                                      -
                                      fe.get_nonzero_components(cell_i).begin());
                       }
-                    
+
                     if (component_mask[component] == true)
                       boundary_values[face_dofs[i]] = dof_values_system[i](component);
-                  } 
+                  }
              }
            else
                                               // fe has only one component,
@@ -1660,9 +1674,9 @@ interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>
                 dof_values_scalar.resize (fe.dofs_per_face);
                function_map.find(boundary_component)->second
                   ->value_list (dof_locations, dof_values_scalar, 0);
-               
+
                                                 // enter into list
-               
+
                for (unsigned int i=0; i<face_dofs.size(); ++i)
                  boundary_values[face_dofs[i]] = dof_values_scalar[i];
              }
@@ -1688,7 +1702,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
 }
 
 #endif
-  
+
 
 template <class DH>
 void
@@ -1699,7 +1713,7 @@ VectorTools::interpolate_boundary_values (const DH                 &dof,
                                          const std::vector<bool>       &component_mask)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping, 
+  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
                              dof, boundary_component,
                              boundary_function, boundary_values, component_mask);
 }
@@ -1714,7 +1728,7 @@ VectorTools::interpolate_boundary_values (const DH                 &dof,
                                          const std::vector<bool>       &component_mask)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping, 
+  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
                              dof, function_map,
                              boundary_values, component_mask);
 }
@@ -1754,7 +1768,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
          ExcMessage ("The number of components in the mask has to be either "
                      "zero or equal to the number of components in the finite "
                      "element."));
-    
+
                                   // check whether boundary values at
                                   // the left or right boundary of
                                   // the line are
@@ -1765,7 +1779,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
                                   // right.
   const unsigned int direction = boundary_component;
   Assert (direction < 2, ExcInvalidBoundaryIndicator());
-  
+
                                   // first find the outermost active
                                   // cell by first traversing the coarse
                                   // grid to its end and then going
@@ -1773,7 +1787,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
   typename DH::cell_iterator outermost_cell = dof.begin(0);
   while (outermost_cell->neighbor(direction).state() == IteratorState::valid)
     outermost_cell = outermost_cell->neighbor(direction);
-  
+
   while (outermost_cell->has_children())
     outermost_cell = outermost_cell->child(direction);
 
@@ -1810,7 +1824,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
   else
     boundary_function.vector_value (outermost_cell->vertex(direction),
                                    function_values);
-  
+
   for (unsigned int i=0; i<fe.dofs_per_vertex; ++i)
     if (component_mask[fe.face_system_to_component_index(i).first])
       {
@@ -1818,7 +1832,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
                                   // entries in the line here?
        const unsigned int row = outermost_cell->vertex_dof_index(direction,i);
        constraints.add_line (row);
-       constraints.set_inhomogeneity (row, 
+       constraints.set_inhomogeneity (row,
                    function_values(fe.face_system_to_component_index(i).first));
       }
 }
@@ -1830,7 +1844,7 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
 // Implementation for 1D
 template <class DH>
 void
-VectorTools::interpolate_boundary_values 
+VectorTools::interpolate_boundary_values
   (const Mapping<DH::dimension, DH::space_dimension>     &mapping,
    const DH                                              &dof,
    const typename FunctionMap<DH::space_dimension>::type &function_map,
@@ -1853,7 +1867,7 @@ VectorTools::interpolate_boundary_values
 
 template <class DH>
 void
-VectorTools::interpolate_boundary_values 
+VectorTools::interpolate_boundary_values
  (const Mapping<DH::dimension, DH::space_dimension>     &mapping,
   const DH                                              &dof,
   const typename FunctionMap<DH::space_dimension>::type &function_map,
@@ -1874,7 +1888,7 @@ VectorTools::interpolate_boundary_values
                                   // immediately
   if (function_map.size() == 0)
     return;
-  
+
   Assert (function_map.find(255) == function_map.end(),
          ExcInvalidBoundaryIndicator());
 
@@ -1922,7 +1936,7 @@ VectorTools::interpolate_boundary_values
   for (unsigned int f=0; f<finite_elements.size(); ++f)
     {
       const FiniteElement<dim> &fe = finite_elements[f];
-      
+
                                       // generate a quadrature rule on the
                                       // face from the unit support
                                       // points. this will be used to
@@ -1973,9 +1987,9 @@ VectorTools::interpolate_boundary_values
              if (component_mask[fe.face_system_to_component_index(i).first]
                  == true)
                unit_support_points[i] = fe.unit_face_support_point(i);
-       
+
          q_collection.push_back (Quadrature<dim-1>(unit_support_points));
-        }    
+        }
     }
                                   // now that we have a q_collection object
                                   // with all the right quadrature points,
@@ -1985,7 +1999,7 @@ VectorTools::interpolate_boundary_values
   hp::MappingCollection<dim> mapping_collection (mapping);
   hp::FEFaceValues<dim> x_fe_values (mapping_collection, finite_elements, q_collection,
                                     update_quadrature_points);
-  
+
   typename DH::active_cell_iterator cell = dof.begin_active(),
                                    endc = dof.end();
   for (; cell!=endc; ++cell)
@@ -2017,10 +2031,10 @@ VectorTools::interpolate_boundary_values
                                    "values that correspond to primitive (scalar) base "
                                    "elements"));
          }
-       
+
        typename DH::face_iterator face = cell->face(face_no);
        const unsigned char boundary_component = face->boundary_indicator();
-       if (function_map.find(boundary_component) != function_map.end()) 
+       if (function_map.find(boundary_component) != function_map.end())
          {
                                             // face is of the right
                                             // component
@@ -2034,7 +2048,7 @@ VectorTools::interpolate_boundary_values
            face->get_dof_indices (face_dofs, cell->active_fe_index());
            const std::vector<Point<DH::space_dimension> > &dof_locations
               = fe_values.get_quadrature_points ();
-           
+
            if (fe_is_system)
              {
                                                  // resize array. avoid
@@ -2046,10 +2060,10 @@ VectorTools::interpolate_boundary_values
                                             Vector<double>(fe.n_components()));
                 else
                   dof_values_system.resize (fe.dofs_per_face);
-                
+
                function_map.find(boundary_component)->second
                   ->vector_value_list (dof_locations, dof_values_system);
-               
+
                                                 // enter those dofs into
                                                 // the list that match the
                                                 // component
@@ -2136,7 +2150,7 @@ VectorTools::interpolate_boundary_values
                                      -
                                      fe.get_nonzero_components(cell_i).begin());
                       }
-                    
+
                     if (component_mask[component] == true)
                      {
                                   // TODO: check whether we should clear
@@ -2145,7 +2159,7 @@ VectorTools::interpolate_boundary_values
                        constraints.set_inhomogeneity (face_dofs[i],
                                                       dof_values_system[i](component));
                      }
-                  } 
+                  }
              }
            else
                                               // fe has only one component,
@@ -2157,9 +2171,9 @@ VectorTools::interpolate_boundary_values
                 dof_values_scalar.resize (fe.dofs_per_face);
                function_map.find(boundary_component)->second
                   ->value_list (dof_locations, dof_values_scalar, 0);
-               
+
                                                 // enter into list
-               
+
                for (unsigned int i=0; i<face_dofs.size(); ++i)
                                   // TODO: check whether we should clear
                                   // the current line first...
@@ -2177,7 +2191,7 @@ VectorTools::interpolate_boundary_values
 
 template <class DH>
 void
-VectorTools::interpolate_boundary_values 
+VectorTools::interpolate_boundary_values
   (const Mapping<DH::dimension, DH::space_dimension> &mapping,
    const DH                                          &dof,
    const unsigned char                                boundary_component,
@@ -2192,11 +2206,11 @@ VectorTools::interpolate_boundary_values
 }
 
 #endif
-  
+
 
 template <class DH>
 void
-VectorTools::interpolate_boundary_values 
+VectorTools::interpolate_boundary_values
   (const DH                            &dof,
    const unsigned char                  boundary_component,
    const Function<DH::space_dimension> &boundary_function,
@@ -2204,7 +2218,7 @@ VectorTools::interpolate_boundary_values
    const std::vector<bool>             &component_mask)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping, 
+  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
                              dof, boundary_component,
                              boundary_function, constraints, component_mask);
 }
@@ -2213,14 +2227,14 @@ VectorTools::interpolate_boundary_values
 
 template <class DH>
 void
-VectorTools::interpolate_boundary_values 
+VectorTools::interpolate_boundary_values
   (const DH                                              &dof,
    const typename FunctionMap<DH::space_dimension>::type &function_map,
    ConstraintMatrix                                      &constraints,
    const std::vector<bool>                               &component_mask)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping, 
+  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
                              dof, function_map,
                              constraints, component_mask);
 }
@@ -2279,16 +2293,16 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
     }
   else
     AssertDimension (dof.get_fe().n_components(), component_mapping.size());
-  
+
   std::vector<unsigned int> dof_to_boundary_mapping;
   std::set<unsigned char> selected_boundary_components;
   for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
        i!=boundary_functions.end(); ++i)
     selected_boundary_components.insert (i->first);
-  
+
   DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components,
                                         dof_to_boundary_mapping);
-  
+
                                   // Done if no degrees of freedom on
                                   // the boundary
   if (dof.n_boundary_dofs (boundary_functions) == 0)
@@ -2318,7 +2332,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
   if (dim>=3)
     {
 #ifdef DEBUG
-// Assert that there are no hanging nodes at the boundary      
+// Assert that there are no hanging nodes at the boundary
       int level = -1;
       for (typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
           cell != dof.end(); ++cell)
@@ -2329,7 +2343,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
                if (level == -1)
                  level = cell->level();
                else
-                 {                       
+                 {
                    Assert (level == cell->level(), ExcNotImplemented());
                  }
              }
@@ -2337,14 +2351,14 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
 #endif
     }
   sparsity.compress();
-  
+
 
                                   // make mass matrix and right hand side
   SparseMatrix<double> mass_matrix(sparsity);
   Vector<double>       rhs(sparsity.n_rows());
 
 
-  MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, 
+  MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
                                              mass_matrix, boundary_functions,
                                              rhs, dof_to_boundary_mapping, (const Function<spacedim>*) 0,
                                              component_mapping);
@@ -2357,7 +2371,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
                                   // eliminate them here.
 
 //TODO: Maybe we should figure out if the element really needs this
-  
+
   FilteredMatrix<Vector<double> > filtered_mass_matrix(mass_matrix, true);
   FilteredMatrix<Vector<double> > filtered_precondition;
   std::vector<bool> excluded_dofs(mass_matrix.m(), false);
@@ -2366,7 +2380,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
   for (unsigned int i=0;i<mass_matrix.m();++i)
     if (mass_matrix.diag_element(i) > max_element)
       max_element = mass_matrix.diag_element(i);
-  
+
   for (unsigned int i=0;i<mass_matrix.m();++i)
     if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
       {
@@ -2375,7 +2389,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
        mass_matrix.diag_element(i) = 1.;
        excluded_dofs[i] = true;
       }
-  
+
   Vector<double> boundary_projection (rhs.size());
 
                                   // Allow for a maximum of 5*n
@@ -2394,7 +2408,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
                                   // solve
   cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
   filtered_precondition.apply_constraints(boundary_projection, true);
-  filtered_precondition.clear();  
+  filtered_precondition.clear();
                                   // fill in boundary values
   for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
     if (dof_to_boundary_mapping[i] != DoFHandler<dim,spacedim>::invalid_dof_index
@@ -2477,16 +2491,16 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
     }
   else
     AssertDimension (dof.get_fe().n_components(), component_mapping.size());
-  
+
   std::vector<unsigned int> dof_to_boundary_mapping;
   std::set<unsigned char> selected_boundary_components;
   for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
        i!=boundary_functions.end(); ++i)
     selected_boundary_components.insert (i->first);
-  
+
   DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components,
                                         dof_to_boundary_mapping);
-  
+
                                   // Done if no degrees of freedom on
                                   // the boundary
   if (dof.n_boundary_dofs (boundary_functions) == 0)
@@ -2516,7 +2530,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
   if (dim>=3)
     {
 #ifdef DEBUG
-// Assert that there are no hanging nodes at the boundary      
+// Assert that there are no hanging nodes at the boundary
       int level = -1;
       for (typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
           cell != dof.end(); ++cell)
@@ -2527,7 +2541,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
                if (level == -1)
                  level = cell->level();
                else
-                 {                       
+                 {
                    Assert (level == cell->level(), ExcNotImplemented());
                  }
              }
@@ -2535,14 +2549,14 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
 #endif
     }
   sparsity.compress();
-  
+
 
                                   // make mass matrix and right hand side
   SparseMatrix<double> mass_matrix(sparsity);
   Vector<double>       rhs(sparsity.n_rows());
 
 
-  MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, 
+  MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
                                              mass_matrix, boundary_functions,
                                              rhs, dof_to_boundary_mapping, (const Function<spacedim>*) 0,
                                              component_mapping);
@@ -2555,7 +2569,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
                                   // eliminate them here.
 
 //TODO: Maybe we should figure out if the element really needs this
-  
+
   FilteredMatrix<Vector<double> > filtered_mass_matrix(mass_matrix, true);
   FilteredMatrix<Vector<double> > filtered_precondition;
   std::vector<bool> excluded_dofs(mass_matrix.m(), false);
@@ -2564,7 +2578,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
   for (unsigned int i=0;i<mass_matrix.m();++i)
     if (mass_matrix.diag_element(i) > max_element)
       max_element = mass_matrix.diag_element(i);
-  
+
   for (unsigned int i=0;i<mass_matrix.m();++i)
     if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
       {
@@ -2573,7 +2587,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
        mass_matrix.diag_element(i) = 1.;
        excluded_dofs[i] = true;
       }
-  
+
   Vector<double> boundary_projection (rhs.size());
 
                                   // Allow for a maximum of 5*n
@@ -2592,7 +2606,7 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
                                   // solve
   cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
   filtered_precondition.apply_constraints(boundary_projection, true);
-  filtered_precondition.clear();  
+  filtered_precondition.clear();
                                   // fill in boundary values
   for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
     if (dof_to_boundary_mapping[i] != DoFHandler<dim,spacedim>::invalid_dof_index
@@ -2650,8 +2664,8 @@ namespace internal
            for (unsigned int i=0; i<dim; ++i)
              dof_indices[i] = numbers::invalid_unsigned_int;
          }
-       
-       
+
+
        bool operator < (const VectorDoFTuple<dim> &other) const
          {
            for (unsigned int i=0; i<dim; ++i)
@@ -2696,7 +2710,7 @@ namespace internal
     add_constraint (const VectorDoFTuple<dim> &dof_indices,
                    const Tensor<1,dim>       &constraining_vector,
                    ConstraintMatrix          &constraints)
-    {      
+    {
                                       // choose the DoF that has the
                                       // largest component in the
                                       // constraining_vector as the
@@ -2748,7 +2762,7 @@ namespace internal
                (std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[2])))
              {
                constraints.add_line (dof_indices.dof_indices[0]);
-               
+
                if (std::fabs (constraining_vector[1]/constraining_vector[0])
                    > std::numeric_limits<double>::epsilon())
                  constraints.add_entry (dof_indices.dof_indices[0],
@@ -2796,14 +2810,14 @@ namespace internal
                                           dof_indices.dof_indices[1],
                                           -constraining_vector[1]/constraining_vector[2]);
                }
-           
+
            break;
          }
 
          default:
                Assert (false, ExcNotImplemented());
        }
-      
+
     }
 
 
@@ -2819,7 +2833,7 @@ namespace internal
     void
     compute_orthonormal_vectors (const Tensor<1,dim> &vector,
                                 Tensor<1,dim> (&orthonormals)[dim-1])
-    {      
+    {
       switch (dim)
        {
          case 3:
@@ -2906,9 +2920,9 @@ namespace internal
 
          default:
                Assert (false, ExcNotImplemented());
-       }      
+       }
     }
-    
+
   }
 }
 
@@ -2926,7 +2940,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
          ExcMessage ("This function is not useful in 1d because it amounts "
                      "to imposing Dirichlet values on the vector-valued "
                      "quantity."));
-         
+
   const FiniteElement<dim> &fe = dof_handler.get_fe();
 
   std::vector<unsigned int> face_dofs (fe.dofs_per_face);
@@ -2949,7 +2963,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                                   // of the same cell. consequently, we also
                                   // have to store which cell a normal vector
                                   // was computed on
-  typedef 
+  typedef
     std::multimap<internal::VectorTools::VectorDoFTuple<dim>,
     std::pair<Tensor<1,dim>, typename DH<dim,spacedim>::active_cell_iterator> >
     DoFToNormalsMap;
@@ -2976,7 +2990,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                                           // vectors at the locations
                                           // where they are defined:
          const std::vector<Point<dim-1> > &
-           unit_support_points = fe.get_unit_face_support_points();  
+           unit_support_points = fe.get_unit_face_support_points();
          Quadrature<dim-1> aux_quad (unit_support_points);
          FEFaceValues<dim> fe_values (mapping, fe, aux_quad,
                                       update_normal_vectors);
@@ -2994,7 +3008,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                                                 // components of vector
                internal::VectorTools::VectorDoFTuple<dim> vector_dofs;
                vector_dofs.dof_indices[0] = face_dofs[i];
-               
+
                for (unsigned int k=0; k<fe.dofs_per_face; ++k)
                  if ((k != i)
                      &&
@@ -3008,7 +3022,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                    vector_dofs.dof_indices[fe.face_system_to_component_index(k).first -
                                            first_vector_component]
                      = face_dofs[k];
-         
+
                for (unsigned int d=0; d<dim; ++d)
                  Assert (vector_dofs.dof_indices[d] < dof_handler.n_dofs(),
                          ExcInternalError());
@@ -3036,7 +3050,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                                   // dofs
   typename DoFToNormalsMap::const_iterator
     p = dof_to_normals_map.begin();
-  
+
   while (p != dof_to_normals_map.end())
     {
                                       // first find the range of entries in
@@ -3128,7 +3142,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                                           // average over all normal vectors
          case 1:
          {
-           
+
                                             // compute the average
                                             // normal vector from all
                                             // the ones that have the
@@ -3145,12 +3159,12 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                   x = cell_to_normals_map.begin();
                 x != cell_to_normals_map.end(); ++x)
              normal += x->second.first;
-           normal /= normal.norm();        
-           
+           normal /= normal.norm();
+
                                             // then construct constraints
                                             // from this:
            const internal::VectorTools::VectorDoFTuple<dim> &
-             dof_indices = same_dof_range[0]->first;       
+             dof_indices = same_dof_range[0]->first;
            internal::VectorTools::add_constraint (dof_indices, normal,
                                                   constraints);
 
@@ -3215,11 +3229,11 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
              for (unsigned int i=0; i<dim; ++i, ++x)
                for (unsigned int j=0; j<dim; ++j)
                  t[i][j] = x->second.first[j];
-             
+
              Assert (std::fabs(determinant (t)) > 1e-3,
                      ExcMessage("Found a set of normal vectors that are nearly collinear."));
            }
-           
+
                                             // so all components of
                                             // this vector dof are
                                             // constrained. enter
@@ -3233,7 +3247,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
 
            break;
          }
-         
+
 
                                           // this is the case of an
                                           // edge contribution in 3d,
@@ -3270,7 +3284,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                 q != same_dof_range[1]; ++q)
              cell_contributions[q->second.second].push_back (q->second.first);
            Assert (cell_contributions.size() >= 1, ExcInternalError());
-           
+
                                             // now for each cell that
                                             // has contributed
                                             // determine the number
@@ -3369,7 +3383,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
 
                Assert (std::fabs (tangent.norm()-1) < 1e-12,
                        ExcInternalError());
-               
+
                tangential_vectors.push_back (tangent);
              }
 
@@ -3392,7 +3406,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH<dim,spacedim>
                if (*t * first_tangent < 0)
                  *t *= -1;
            }
-           
+
                                             // now compute the
                                             // average tangent and
                                             // normalize it
@@ -3457,7 +3471,7 @@ namespace internal
                                        // necessary, so have a read-write
                                        // version of it:
       double exponent = exponent_1;
-  
+
       const unsigned int        n_components = dof.get_fe().n_components();
       const bool                fe_is_system = (n_components != 1);
 
@@ -3468,7 +3482,7 @@ namespace internal
         }
 
       difference.reinit (dof.get_tria().n_active_cells());
-  
+
       switch (norm)
         {
           case dealii::VectorTools::L2_norm:
@@ -3482,7 +3496,7 @@ namespace internal
           default:
                 break;
         }
-  
+
       UpdateFlags update_flags = UpdateFlags (update_quadrature_points  |
                                               update_JxW_values);
       switch (norm)
@@ -3500,13 +3514,13 @@ namespace internal
           default:
                 update_flags |= UpdateFlags (update_values);
                 break;
-        }  
+        }
 
       dealii::hp::FECollection<dim,spacedim> fe_collection (dof.get_fe());
       dealii::hp::FEValues<dim,spacedim> x_fe_values(mapping, fe_collection, q, update_flags);
 
       const unsigned int max_n_q_points = q.max_n_quadrature_points ();
-      
+
       std::vector< dealii::Vector<double> >
         function_values (max_n_q_points, dealii::Vector<double>(n_components));
       std::vector<std::vector<Tensor<1,spacedim> > >
@@ -3529,7 +3543,7 @@ namespace internal
                                        // scalar functions
       std::vector<double>         tmp_values (max_n_q_points);
       std::vector<Tensor<1,spacedim> > tmp_gradients (max_n_q_points);
-  
+
                                        // loop over all cells
       typename DH::active_cell_iterator cell = dof.begin_active(),
                                         endc = dof.end();
@@ -3541,7 +3555,7 @@ namespace internal
 
           const dealii::FEValues<dim, spacedim> &fe_values  = x_fe_values.get_present_fe_values ();
           const unsigned int   n_q_points = fe_values.n_quadrature_points;
-          
+
                                            // resize all out scratch
                                            // arrays to the number of
                                            // quadrature points we use
@@ -3582,8 +3596,8 @@ namespace internal
               for (unsigned int k=0;k<n_q_points;++k)
                 weight_vectors[k] = 1.;
             }
-      
-      
+
+
           if (update_flags & update_values)
             {
                                                // first compute the exact solution
@@ -3603,7 +3617,7 @@ namespace internal
                   for (unsigned int i=0; i<n_q_points; ++i)
                     psi_values[i](0) = tmp_values[i];
                 }
-         
+
                                                // then subtract finite element
                                                // fe_function
               fe_values.get_function_values (fe_function, function_values);
@@ -3629,8 +3643,8 @@ namespace internal
                                                 tmp_gradients);
                   for (unsigned int i=0; i<n_q_points; ++i)
                     psi_grads[i][0] = tmp_gradients[i];
-                }      
-         
+                }
+
                                                // then subtract finite element
                                                // function_grads
               fe_values.get_function_grads (fe_function, function_grads);
@@ -3638,7 +3652,7 @@ namespace internal
                 for (unsigned int q=0; q<n_q_points; ++q)
                   psi_grads[q][k] -= function_grads[q][k];
             }
-      
+
           switch (norm)
             {
               case dealii::VectorTools::mean:
@@ -3666,7 +3680,7 @@ namespace internal
                         psi_scalar[q] += std::pow(psi_values[q](k)*psi_values[q](k),
                                                   exponent/2.)
                                          * weight_vectors[q](k);
-           
+
                                                      // Integrate
                     diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
                                                fe_values.get_JxW_values().begin(),
@@ -3730,7 +3744,7 @@ namespace internal
                         psi_scalar[q] += std::pow(psi_grads[q][k] * psi_grads[q][k],
                                                   exponent/2.)
                                          * weight_vectors[q](k);
-           
+
                     diff += std::inner_product (psi_scalar.begin(), psi_scalar.end(),
                                                 fe_values.get_JxW_values().begin(),
                                                 0.0);
@@ -3763,7 +3777,7 @@ namespace internal
                           for (unsigned int d=0;d<dim;++d)
                             t = std::max(t,std::fabs(psi_grads[q][k][d])
                                          * weight_vectors[q](k));
-                 
+
                           psi_scalar[q] = std::max(psi_scalar[q],t);
                         }
 
@@ -3907,7 +3921,7 @@ VectorTools::point_difference (const Mapping<dim, spacedim>    &mapping,
 
   Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
          ExcInternalError());
-  
+
   const Quadrature<dim>
     quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
   FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
@@ -3922,7 +3936,7 @@ VectorTools::point_difference (const Mapping<dim, spacedim>    &mapping,
     difference(0) = exact_function.value(point);
   else
     exact_function.vector_value(point, difference);
-    
+
   for (unsigned int i=0; i<difference.size(); ++i)
     difference(i) -= u_value[0](i);
 }
@@ -4016,7 +4030,7 @@ VectorTools::point_value (const Mapping<dim, spacedim>    &mapping,
 
   Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
          ExcInternalError());
-  
+
   const Quadrature<dim>
     quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
   FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
@@ -4043,7 +4057,7 @@ VectorTools::compute_mean_value (const Mapping<dim, spacedim>    &mapping,
          ExcDimensionMismatch (v.size(), dof.n_dofs()));
   Assert (component < dof.get_fe().n_components(),
          ExcIndexRange(component, 0, dof.get_fe().n_components()));
-  
+
   FEValues<dim> fe(mapping, dof.get_fe(), quadrature,
                   UpdateFlags(update_JxW_values
                               | update_values));
@@ -4051,7 +4065,7 @@ VectorTools::compute_mean_value (const Mapping<dim, spacedim>    &mapping,
   typename DoFHandler<dim,spacedim>::active_cell_iterator c;
   std::vector<Vector<double> > values(quadrature.size(),
                                      Vector<double> (dof.get_fe().n_components()));
-  
+
   double mean = 0.;
   double area = 0.;
                                   // Compute mean value
@@ -4065,7 +4079,7 @@ VectorTools::compute_mean_value (const Mapping<dim, spacedim>    &mapping,
          area += fe.JxW(k);
        };
     };
-  
+
   return (mean/area);
 }
 
index 58d8142dc2573ef2e3bd60ea34670f8e822f2a61..5b724676d02c614701da3009bef1c1d14b0ee5f3 100644 (file)
@@ -2018,6 +2018,7 @@ void MatrixCreator::create_mass_matrix<deal_II_dimension>
  const Function<deal_II_dimension> * const coefficient);
 
 
+#if deal_II_dimension != 1
 
 template
 void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
@@ -2055,7 +2056,6 @@ MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
  const Function<deal_II_dimension> * const,
  std::vector<unsigned int>);
 
-
 template
 void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
 (const hp::DoFHandler<deal_II_dimension>&,
@@ -2067,6 +2067,7 @@ void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
  const Function<deal_II_dimension> * const,
  std::vector<unsigned int>);
 
+#endif
 
 
 // non-hp versions of create_laplace_matrix
@@ -2201,6 +2202,8 @@ void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_dimension
  Vector<double>           &rhs_vector,
  const Function<deal_II_dimension+1> * const coefficient);
 
+#if deal_II_dimension == 2
+
 template
 void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_dimension+1>
 (const DoFHandler<deal_II_dimension,deal_II_dimension+1>     &dof,
@@ -2224,7 +2227,7 @@ void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_dimens
  const Function<deal_II_dimension+1> * const a,
  std::vector<unsigned int>);
 
+#endif 
 
 // #if deal_II_dimension != 1
 // template

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.