From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 30 Aug 2005 09:48:05 +0000 (+0000)
Subject: move instatiations to vectors.templates.h; thanks to Dima Sorkin
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=846311e3b77a953d15f3035649fc735c8709d899;p=dealii-svn.git

move instatiations to vectors.templates.h; thanks to Dima Sorkin


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

diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h
index 22447c2012..2b37acf61a 100644
--- a/deal.II/deal.II/include/numerics/vectors.h
+++ b/deal.II/deal.II/include/numerics/vectors.h
@@ -292,6 +292,7 @@ class ConstraintMatrix;
  * if access to an object describing the exact form of the boundary is needed, the
  * pointer stored within the triangulation object is accessed.
  *
+ * @ref Instantiations: some (Vector<float>, Vector<double>, BlockVector<float>, BlockVector<double>, see also individual functions)
  * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, 1998, 1999, 2000, 2001
  */
 class VectorTools
@@ -432,6 +433,15 @@ class VectorTools
 				      * space @p dof_2 afterwards, to
 				      * make the result continuous
 				      * again.
+				      *
+				      * @ref Instantiations: some (the
+				      * standard vector types
+				      * mentioned in the class
+				      * documentation are
+				      * instantiated, but only the
+				      * same vector for InVector and
+				      * OutVector. Other combinations
+				      * must be instantiated by hand.)
 				      */
     template <int dim, class InVector, class OutVector>
     static void interpolate (const DoFHandler<dim>    &dof_1,
@@ -830,6 +840,12 @@ class VectorTools
 				      *
 				      * See the general documentation of this
 				      * class for more information.
+				      *
+				      * @ref Instantiations: some
+				      * (InVectors as in the
+				      * documentation of the class,
+				      * OutVector only Vector<double>
+				      * and Vector<float>)
 				      */
     template <int dim, class InVector, class OutVector>
     static void integrate_difference (const Mapping<dim>    &mapping,
diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h
new file mode 100644
index 0000000000..4f07b685c2
--- /dev/null
+++ b/deal.II/deal.II/include/numerics/vectors.templates.h
@@ -0,0 +1,1648 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef _deal2__vectors_templates_h
+#define _deal2__vectors_templates_h
+
+#include <base/function.h>
+#include <base/quadrature.h>
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/precondition.h>
+#include <lac/solver_cg.h>
+#include <lac/vector_memory.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_tools.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+
+#include <numeric>
+#include <algorithm>
+#include <vector>
+#include <cmath>
+
+template <int dim>
+inline double sqr_point (const Tensor<1,dim> &p)
+{
+  return p * p;
+}
+
+
+template <int dim, class VECTOR>
+void VectorTools::interpolate (const Mapping<dim>    &mapping,
+			       const DoFHandler<dim> &dof,
+			       const Function<dim>   &function,
+			       VECTOR                &vec)
+{
+  Assert (dof.get_fe().n_components() == function.n_components,
+	  ExcComponentMismatch());
+  
+  const FiniteElement<dim> &fe           = dof.get_fe();
+  const unsigned int        n_components = fe.n_components();
+  const bool                fe_is_system = (n_components != 1);
+  
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+						 endc = dof.end();
+
+				   // For FESystems many of the
+				   // unit_support_points will
+				   // appear multiply, as a point
+				   // may be unit_support_point
+				   // for several of the components
+				   // of the system.
+				   // The following is rather
+				   // complicated as it is
+				   // avoided to evaluate
+				   // the vectorfunction multiply at
+				   // the same point on a cell.
+  const std::vector<Point<dim> > &
+    unit_support_points = fe.get_unit_support_points();
+  Assert (unit_support_points.size() != 0,
+	  ExcNonInterpolatingFE());
+
+				   // Find the support points 
+				   // on a cell that
+				   // are multiply mentioned in 
+				   // @p{unit_support_points}.
+				   // Mark the first representative
+				   // of each multiply mentioned
+				   // support point by appending its
+				   // dof index to @p{dofs_of_rep_points}.
+				   // Each multiple point gets to know
+				   // the dof index of its representative
+				   // point by the @p{dof_to_rep_dof_table}.
+
+				   // the following vector collects all dofs i,
+				   // 0<=i<fe.dofs_per_cell, for that
+				   // 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
+				   // is called rep index.
+  std::vector<unsigned int> dofs_of_rep_points;
+				   // the following table converts a dof i
+				   // to the rep index.
+  std::vector<unsigned int> dof_to_rep_index_table;
+  unsigned int n_rep_points=0;
+  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+    {
+      bool representative=true;
+				       // the following loop is looped
+				       // the other way round to get
+				       // the minimal effort of
+				       // O(fe.dofs_per_cell) for multiple
+				       // support points that are placed
+				       // one after the other.
+      for (unsigned int j=dofs_of_rep_points.size(); j>0; --j)
+	if (unit_support_points[i] 
+	    == unit_support_points[dofs_of_rep_points[j-1]])
+	  {
+	    dof_to_rep_index_table.push_back(j-1);
+	    representative=false;
+	    break;
+	  }
+      
+      if (representative)
+	{
+					   // rep_index=dofs_of_rep_points.size()
+	  dof_to_rep_index_table.push_back(dofs_of_rep_points.size());
+					   // dofs_of_rep_points[rep_index]=i
+	  dofs_of_rep_points.push_back(i);
+	  ++n_rep_points;
+	}
+    }
+  Assert(dofs_of_rep_points.size()==n_rep_points, ExcInternalError());
+  Assert(dof_to_rep_index_table.size()==fe.dofs_per_cell, ExcInternalError());
+
+  std::vector<unsigned int> dofs_on_cell (fe.dofs_per_cell);
+  std::vector<Point<dim> >  rep_points (n_rep_points);
+
+				   // get space for the values of the
+				   // function at the rep support points.
+				   //
+				   // have two versions, one for system fe
+				   // and one for scalar ones, to take the
+				   // more efficient one respectively
+  std::vector<double>          function_values_scalar (n_rep_points);
+  std::vector<Vector<double> > function_values_system (n_rep_points,
+						  Vector<double>(fe.n_components()));
+
+				   // Make a quadrature rule from support points
+				   // to feed it into FEValues
+  Quadrature<dim> support_quadrature(unit_support_points);
+
+				   // Transformed support points are computed by
+				   // FEValues
+  FEValues<dim> fe_values (mapping, fe, support_quadrature, update_q_points);
+  
+  for (; cell!=endc; ++cell)
+    {
+				       // for each cell:
+				       // get location of finite element
+				       // support_points
+      fe_values.reinit(cell);
+      const std::vector<Point<dim> >& support_points =
+	fe_values.get_quadrature_points();
+      
+				       // pick out the representative
+				       // support points
+      for (unsigned int j=0; j<dofs_of_rep_points.size(); ++j)
+	rep_points[j]=support_points[dofs_of_rep_points[j]];
+
+				       // get indices of the dofs on this cell
+      cell->get_dof_indices (dofs_on_cell);
+
+
+      if (fe_is_system)
+	{
+					   // get function values at
+					   // these points. Here: get
+					   // all components
+	  function.vector_value_list (rep_points, function_values_system);
+					   // distribute the function
+					   // values to the global
+					   // vector
+	  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+	    {
+	      const unsigned int component
+		= fe.system_to_component_index(i).first;
+	      const unsigned int rep_dof=dof_to_rep_index_table[i];
+	      vec(dofs_on_cell[i])
+		= function_values_system[rep_dof](component);
+	    };
+	}
+      
+      else
+	{
+					   // get first component only,
+					   // which is the only component
+					   // in the function anyway
+	  function.value_list (rep_points, function_values_scalar, 0);
+					   // distribute the function
+					   // values to the global
+					   // vector
+	  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+	    vec(dofs_on_cell[i]) 
+	      = function_values_scalar[dof_to_rep_index_table[i]];
+	};
+    }
+}
+
+
+template <int dim, class VECTOR>
+void VectorTools::interpolate (const DoFHandler<dim> &dof,
+			       const Function<dim>   &function,
+			       VECTOR                &vec)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  interpolate(mapping, dof, function, vec);
+}
+
+
+
+
+template <int dim, class InVector, class OutVector>
+void
+VectorTools::interpolate (const DoFHandler<dim>           &dof_1,
+			  const DoFHandler<dim>           &dof_2,
+			  const FullMatrix<double>        &transfer,
+			  const InVector                  &data_1,
+			  OutVector                       &data_2)
+{
+  Vector<double> cell_data_1(dof_1.get_fe().dofs_per_cell);
+  Vector<double> cell_data_2(dof_2.get_fe().dofs_per_cell);
+
+  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>::active_cell_iterator h = dof_1.begin_active();
+  typename DoFHandler<dim>::active_cell_iterator l = dof_2.begin_active();
+  const typename DoFHandler<dim>::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) 
+      {
+	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());	
+	++touch_count[local_dof_indices[j]];
+      };
+  };
+
+				   // compute the mean value of the
+				   // sum which we have placed in each
+				   // entry of the output vector
+  for (unsigned int i=0; i<dof_2.n_dofs(); ++i)
+    {
+      Assert (touch_count[i] != 0,
+	      ExcInternalError());
+      
+      data_2(i) /= touch_count[i];
+    };
+}
+
+
+#if deal_II_dimension == 1
+
+void VectorTools::project (const Mapping<1>       &,
+			   const DoFHandler<1>    &,
+			   const ConstraintMatrix &,
+			   const Quadrature<1>    &,
+			   const Function<1>      &,
+			   Vector<double>         &,
+			   const bool              ,
+			   const Quadrature<0>    &,
+			   const bool              )
+{
+				   // this function should easily be implemented
+				   // using the template below. However some
+				   // changes have to be made since faces don't
+				   // exist in 1D. Maybe integrate the creation of
+				   // zero boundary values into the
+				   // project_boundary_values function?
+  Assert (false, ExcNotImplemented());
+}
+
+
+#endif
+
+
+template <int dim>
+void VectorTools::project (const Mapping<dim>       &mapping,
+			   const DoFHandler<dim>    &dof,
+			   const ConstraintMatrix   &constraints,
+			   const Quadrature<dim>    &quadrature,
+			   const Function<dim>      &function,
+			   Vector<double>           &vec,
+			   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,
+	  ExcInvalidFE());
+  
+  const FiniteElement<dim> &fe = dof.get_fe();
+
+				   // make up boundary values
+  std::map<unsigned int,double> boundary_values;
+
+  if (enforce_zero_boundary == true) 
+				     // no need to project boundary
+				     // values, but enforce
+				     // homogeneous boundary values
+				     // anyway
+    {
+				       // loop over all boundary faces
+				       // to get all dof indices of
+				       // dofs on the boundary. note
+				       // that in 3d there are cases
+				       // where a face is not at the
+				       // boundary, yet one of its
+				       // lines is, and we should
+				       // consider the degrees of
+				       // freedom on it as boundary
+				       // nodes. likewise, in 2d and
+				       // 3d there are cases where a
+				       // cell is only at the boundary
+				       // by one vertex. nevertheless,
+				       // since we do not support
+				       // boundaries with dimension
+				       // less or equal to dim-2, each
+				       // such boundary dof is also
+				       // found from some other face
+				       // that is actually wholly on
+				       // the boundary, not only by
+				       // one line or one vertex
+      typename DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+						     endf = dof.end_face();
+      std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
+      for (; face!=endf; ++face)
+	if (face->at_boundary())
+	  {
+	    face->get_dof_indices (face_dof_indices);
+	    for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+					       // enter zero boundary values
+					       // for all boundary nodes
+					       //
+					       // we need not care about
+					       // vector valued elements here,
+					       // since we set all components
+	      boundary_values[face_dof_indices[i]] = 0.;
+	  };
+    }
+  else
+				     // no homogeneous boundary values
+    if (project_to_boundary_first == true)
+				       // boundary projection required
+      {
+					 // set up a list of boundary functions for
+					 // the different boundary parts. We want the
+					 // @p{function} to hold on all parts of the
+					 // boundary
+	typename FunctionMap<dim>::type boundary_functions;
+	for (unsigned char c=0; c<255; ++c)
+	  boundary_functions[c] = &function;
+	project_boundary_values (dof, boundary_functions, q_boundary,
+				 boundary_values);
+      };
+
+
+				   // set up mass matrix and right hand side
+  vec.reinit (dof.n_dofs());
+  SparsityPattern sparsity(dof.n_dofs(),
+			   dof.n_dofs(),
+			   dof.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof, sparsity);
+  constraints.condense (sparsity);
+  
+  SparseMatrix<double> mass_matrix (sparsity);
+  Vector<double> tmp (mass_matrix.n());
+
+  MatrixCreator::create_mass_matrix (mapping, dof, quadrature, mass_matrix);
+  
+  VectorTools::create_right_hand_side (mapping, dof, quadrature, function, tmp);
+
+  constraints.condense (mass_matrix);
+  constraints.condense (tmp);
+  if (boundary_values.size() != 0)
+    MatrixTools::apply_boundary_values (boundary_values,
+					mass_matrix, vec, tmp,
+					true);
+
+  SolverControl           control(1000,1e-16);
+  PrimitiveVectorMemory<> memory;
+  SolverCG<>              cg(control,memory);
+
+  PreconditionSSOR<> prec;
+  prec.initialize(mass_matrix, 1.2);
+				   // solve
+  cg.solve (mass_matrix, vec, tmp, prec);
+  
+				   // distribute solution
+  constraints.distribute (vec);
+}
+
+
+template <int dim>
+void VectorTools::project (const DoFHandler<dim>    &dof,
+			   const ConstraintMatrix   &constraints,
+			   const Quadrature<dim>    &quadrature,
+			   const Function<dim>      &function,
+			   Vector<double>           &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"));
+  static const MappingQ1<dim> mapping;
+  project(mapping, dof, constraints, quadrature, function, vec,
+	  enforce_zero_boundary, q_boundary, project_to_boundary_first);
+}
+
+
+
+
+template <int dim>
+void VectorTools::create_right_hand_side (const Mapping<dim>    &mapping,
+					  const DoFHandler<dim> &dof_handler,
+					  const Quadrature<dim> &quadrature,
+					  const Function<dim>   &rhs_function,
+					  Vector<double>        &rhs_vector)
+{
+  const FiniteElement<dim> &fe  = dof_handler.get_fe();
+  Assert (fe.n_components() == rhs_function.n_components,
+	  ExcComponentMismatch());
+  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_q_points |
+					 update_JxW_values);
+  FEValues<dim> fe_values (mapping, fe, quadrature, update_flags);
+
+  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>::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) 
+	{
+	  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) 
+	      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));
+      
+      // Use the faster code if the FiniteElement is primitive
+      if (fe.is_primitive ())
+	{
+	  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;
+	      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];
+		  }
+	      
+	      cell->get_dof_indices (dofs);
+	      
+	      for (unsigned int i=0; i<dofs_per_cell; ++i)
+		rhs_vector(dofs[i]) += cell_vector(i);
+	    }
+	}
+      else
+	// Otherwise do it the way proposed for vector valued elements
+	{
+	  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;
+	      for (unsigned int point=0; point<n_q_points; ++point)
+		for (unsigned int i=0; i<dofs_per_cell; ++i)
+		  for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+		    if (fe.get_nonzero_components(i)[comp_i])
+		      {
+			cell_vector(i) += rhs_values[point](comp_i) *
+			                  fe_values.shape_value_component(i,point,comp_i) *
+			                  weights[point];
+		      }
+	      
+	      cell->get_dof_indices (dofs);
+	      
+	      for (unsigned int i=0; i<dofs_per_cell; ++i)
+		rhs_vector(dofs[i]) += cell_vector(i);
+	    }
+	}
+    }
+}
+
+
+template <int dim>
+void VectorTools::create_right_hand_side (const DoFHandler<dim>    &dof_handler,
+					  const Quadrature<dim>    &quadrature,
+					  const Function<dim>      &rhs_function,
+					  Vector<double>           &rhs_vector)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  create_right_hand_side(mapping, dof_handler, quadrature,
+			 rhs_function, rhs_vector);
+}
+
+
+
+#if deal_II_dimension != 1
+
+template <int dim>
+void
+VectorTools::create_boundary_right_hand_side (const Mapping<dim>      &mapping,
+					      const DoFHandler<dim>   &dof_handler,
+					      const Quadrature<dim-1> &quadrature,
+					      const Function<dim>     &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,
+	  ExcComponentMismatch());
+  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_q_points |
+					 update_JxW_values);
+  FEFaceValues<dim> fe_values (mapping, fe, quadrature, update_flags);
+
+  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>::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 () &&
+	      (boundary_indicators.find (cell->face(face)->boundary_indicator())
+	       !=
+	       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) 
+		  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 (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())
+	       !=
+	       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];
+		      }
+		}
+	      else
+		{
+		  // And the full featured code, if vector valued FEs are used
+		  for (unsigned int point=0; point<n_q_points; ++point)
+		    for (unsigned int i=0; i<dofs_per_cell; ++i)
+		      for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+			if (fe.get_nonzero_components(i)[comp_i])
+			  {
+			    cell_vector(i) += rhs_values[point](comp_i) *
+			                      fe_values.shape_value_component(i,point,comp_i) *
+			                      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
+
+void
+VectorTools::create_boundary_right_hand_side (const Mapping<1>    &,
+					      const DoFHandler<1> &,
+					      const Quadrature<0> &,
+					      const Function<1>   &,
+					      Vector<double>      &,
+					      const std::set<unsigned char> &)
+{
+  Assert (false, ExcImpossibleInDim(1));
+}
+
+#endif
+
+template <int dim>
+void
+VectorTools::create_boundary_right_hand_side (const DoFHandler<dim>   &dof_handler,
+					      const Quadrature<dim-1> &quadrature,
+					      const Function<dim>     &rhs_function,
+					      Vector<double>          &rhs_vector,
+					      const std::set<unsigned char> &boundary_indicators)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  create_boundary_right_hand_side(mapping, dof_handler, quadrature,
+				  rhs_function, rhs_vector,
+				  boundary_indicators);
+}
+
+
+
+#if deal_II_dimension == 1
+
+void
+VectorTools::interpolate_boundary_values (const Mapping<1>         &,
+					  const DoFHandler<1>      &dof,
+					  const unsigned char       boundary_component,
+					  const Function<1>        &boundary_function,
+					  std::map<unsigned int,double> &boundary_values,
+					  const std::vector<bool>       &component_mask_)
+{
+  Assert (boundary_component != 255,
+	  ExcInvalidBoundaryIndicator());
+
+  const FiniteElement<1> &fe = dof.get_fe();
+  Assert (fe.n_components() == boundary_function.n_components,
+	  ExcComponentMismatch());
+
+				   // set the component mask to either
+				   // the original value or a vector
+				   // of @p{true}s
+  const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
+					  std::vector<bool> (fe.n_components(), true) :
+					  component_mask_);
+  Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
+	  ExcComponentMismatch());
+  
+				   // check whether boundary values at
+				   // the left or right boundary of
+				   // the line are
+				   // requested. @p{direction} denotes
+				   // the neighboring direction in
+				   // which we seek the boundary,
+				   // i.e. 0 is left boundary and 1 is
+				   // 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
+				   // to the children
+  DoFHandler<1>::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);
+
+				   // now set the value of the
+				   // outermost degree of
+				   // freedom. setting also
+				   // creates the entry in the map
+				   // if it did not exist
+				   // beforehand
+				   //
+				   // save some time by requesting
+				   // values only once for each point,
+				   // irrespective of the number of
+				   // components of the function
+  Vector<double> function_values (fe.n_components());
+  if (fe.n_components() == 1)
+    function_values(0)
+      = boundary_function.value (outermost_cell->vertex(direction));
+  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)]
+	= function_values(fe.face_system_to_component_index(i).first);
+}
+
+
+
+void
+VectorTools::interpolate_boundary_values (const Mapping<1>              &mapping,
+					  const DoFHandler<1>           &dof,
+					  const FunctionMap<1>::type    &function_map,
+					  std::map<unsigned int,double> &boundary_values,
+					  const std::vector<bool>       &component_mask)
+{
+  for (FunctionMap<1>::type::const_iterator i=function_map.begin();
+       i!=function_map.end(); ++i)
+    interpolate_boundary_values (mapping, dof, i->first, *i->second,
+				 boundary_values, component_mask);
+}
+
+
+#endif
+
+
+template <int dim>
+void
+VectorTools::
+interpolate_boundary_values (const Mapping<dim>            &mapping,
+                             const DoFHandler<dim>         &dof,
+                             const typename FunctionMap<dim>::type &function_map,
+                             std::map<unsigned int,double> &boundary_values,
+                             const std::vector<bool>       &component_mask_)
+{
+				   // if for whatever reason we were
+				   // passed an empty map, return
+				   // immediately
+  if (function_map.size() == 0)
+    return;
+  
+  Assert (function_map.find(255) == function_map.end(),
+	  ExcInvalidBoundaryIndicator());
+
+  const FiniteElement<dim> &fe           = dof.get_fe();
+  const unsigned int        n_components = fe.n_components();
+  const bool                fe_is_system = (n_components != 1);
+
+  for (typename FunctionMap<dim>::type::const_iterator i=function_map.begin();
+       i!=function_map.end(); ++i)
+    Assert (n_components == i->second->n_components,
+	    ExcInvalidFE());
+
+				   // set the component mask to either
+				   // the original value or a vector
+				   // of @p{true}s
+  const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
+					  std::vector<bool> (fe.n_components(), true) :
+					  component_mask_);
+  Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
+	  ExcComponentMismatch());
+
+				   // field to store the indices
+  std::vector<unsigned int> face_dofs (fe.dofs_per_face,
+				       DoFHandler<dim>::invalid_dof_index);
+  std::vector<Point<dim> >  dof_locations (face_dofs.size(), Point<dim>());
+  
+				   // array to store the values of
+				   // the boundary function at the
+				   // boundary points. have to arrays
+				   // for scalar and vector functions
+				   // to use the more efficient one
+				   // respectively
+  std::vector<double>          dof_values_scalar (fe.dofs_per_face);
+  std::vector<Vector<double> > dof_values_system (fe.dofs_per_face,
+						  Vector<double>(fe.n_components()));
+
+				   // next generate a quadrature rule
+				   // on the face from the unit
+				   // support points. this wil be used
+				   // to obtain the quadrature points
+				   // on the real cell's face
+  std::vector<Point<dim-1> >
+    unit_support_points = fe.get_unit_face_support_points();
+  
+				   // check whether there are support
+				   // points on the face. if not, then
+				   // we should try a more clever
+				   // way. the idea is that a finite
+				   // element may not offer support
+				   // points for all its shape
+				   // functions, but maybe only
+				   // some. if it offers support
+				   // points for the components we are
+				   // interested in in this function,
+				   // then that's fine. if not, the
+				   // function we call in the finite
+				   // element will raise an
+				   // exception. the support points
+				   // for the other shape functions
+				   // are left uninitialized (well,
+				   // initialized by the default
+				   // constructor), since we don't
+				   // need them anyway.
+  if (unit_support_points.size() == 0)
+    {
+      unit_support_points.resize (fe.dofs_per_face);
+      for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+        if (component_mask[fe.face_system_to_component_index(i).first]
+            == true)
+          unit_support_points[i] = fe.unit_face_support_point(i);
+    };
+
+  Quadrature<dim-1> aux_quad (unit_support_points);
+  FEFaceValues<dim> fe_values (mapping, fe, aux_quad, update_q_points);
+
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+						 endc = dof.end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+	 ++face_no)
+      {
+	typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+	const unsigned char boundary_component = face->boundary_indicator();
+	if (function_map.find(boundary_component) != function_map.end()) 
+					   // face is of the right component
+	  {
+					     // get indices, physical location and
+					     // boundary values of dofs on this
+					     // face
+	    face->get_dof_indices (face_dofs);
+	    fe_values.reinit(cell, face_no);
+	    const std::vector<Point<dim> > &dof_locations = fe_values.get_quadrature_points ();
+	    
+	    if (fe_is_system)
+	      {
+		function_map.find(boundary_component)->second
+                  ->vector_value_list (dof_locations, dof_values_system);
+		
+						 // enter those dofs
+						 // into the list that
+						 // match the
+						 // component
+						 // signature. avoid
+						 // the usual
+						 // complication that
+						 // we can't just use
+						 // *_system_to_component_index
+						 // for non-primitive
+						 // FEs
+		for (unsigned int i=0; i<face_dofs.size(); ++i)
+                  {
+                    unsigned int component;
+                    if (fe.is_primitive())
+                      component = fe.face_system_to_component_index(i).first;
+                    else
+                      {
+                                                         // non-primitive
+                                                         // case. make
+                                                         // sure that
+                                                         // this
+                                                         // particular
+                                                         // shape
+                                                         // function
+                                                         // _is_
+                                                         // primitive,
+                                                         // and get at
+                                                         // it's
+                                                         // component. use
+                                                         // usual
+                                                         // trick to
+                                                         // transfer
+                                                         // face dof
+                                                         // index to
+                                                         // cell dof
+                                                         // index
+                        const unsigned int cell_i
+                          = (dim == 1 ?
+                             i
+                             :
+                             (dim == 2 ?
+                              (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+                              :
+                              (dim == 3 ?
+                               (i<4*fe.dofs_per_vertex ?
+                                i
+                                :
+                                (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
+                                 i+4*fe.dofs_per_vertex
+                                 :
+                                 i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+                               :
+                               deal_II_numbers::invalid_unsigned_int)));
+                        Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
+
+                                                         // make sure
+                                                         // that if
+                                                         // this is
+                                                         // not a
+                                                         // primitive
+                                                         // shape function,
+                                                         // then all
+                                                         // the
+                                                         // corresponding
+                                                         // components
+                                                         // in the
+                                                         // mask are
+                                                         // not set
+                        if (!fe.is_primitive(cell_i))
+                          for (unsigned int c=0; c<n_components; ++c)
+                            if (fe.get_nonzero_components(cell_i)[c])
+                              Assert (component_mask[c] == false,
+                                      ExcFENotPrimitive());
+
+                                                         // let's pick
+                                                         // the first
+                                                         // of
+                                                         // possibly
+                                                         // more than
+                                                         // one
+                                                         // non-zero
+                                                         // components. if
+                                                         // shape
+                                                         // function
+                                                         // is
+                                                         // non-primitive,
+                                                         // then we
+                                                         // will
+                                                         // ignore the
+                                                         // result in
+                                                         // the
+                                                         // following
+                                                         // anyway,
+                                                         // otherwise
+                                                         // there's
+                                                         // only one
+                                                         // non-zero
+                                                         // component
+                                                         // which we
+                                                         // will use
+                        component = (std::find (fe.get_nonzero_components(cell_i).begin(),
+                                                fe.get_nonzero_components(cell_i).end(),
+                                                true)
+                                     -
+                                     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,
+					       // so save some computations
+	      {
+						 // get only the one component that
+						 // this function has
+		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];
+	      }
+	  }
+      }
+}
+
+
+
+template <int dim>
+void
+VectorTools::interpolate_boundary_values (const Mapping<dim>            &mapping,
+					  const DoFHandler<dim>         &dof,
+					  const unsigned char            boundary_component,
+					  const Function<dim>           &boundary_function,
+					  std::map<unsigned int,double> &boundary_values,
+					  const std::vector<bool>       &component_mask)
+{
+  typename FunctionMap<dim>::type function_map;
+  function_map[boundary_component] = &boundary_function;
+  interpolate_boundary_values (mapping, dof, function_map, boundary_values,
+			       component_mask);
+}
+
+
+  
+template <int dim>
+void
+VectorTools::interpolate_boundary_values (const DoFHandler<dim>         &dof,
+					  const unsigned char            boundary_component,
+					  const Function<dim>           &boundary_function,
+					  std::map<unsigned int,double> &boundary_values,
+					  const std::vector<bool>       &component_mask)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  interpolate_boundary_values(mapping, dof, boundary_component,
+			      boundary_function, boundary_values, component_mask);
+}
+
+
+
+template <int dim>
+void
+VectorTools::interpolate_boundary_values (const DoFHandler<dim>         &dof,
+					  const typename FunctionMap<dim>::type &function_map,
+					  std::map<unsigned int,double> &boundary_values,
+					  const std::vector<bool>       &component_mask)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  interpolate_boundary_values(mapping, dof, function_map,
+			      boundary_values, component_mask);
+}
+
+
+#if deal_II_dimension == 1
+
+void
+VectorTools::project_boundary_values (const Mapping<1>       &mapping,
+				      const DoFHandler<1>    &dof,
+				      const FunctionMap<1>::type &boundary_functions,
+				      const Quadrature<0>  &,
+				      std::map<unsigned int,double> &boundary_values)
+{
+				   // projection in 1d is equivalent
+				   // to interpolation
+  interpolate_boundary_values (mapping, dof, boundary_functions,
+			       boundary_values, std::vector<bool>());
+}
+
+#endif
+
+
+template <int dim>
+void
+VectorTools::project_boundary_values (const Mapping<dim>       &mapping,
+				      const DoFHandler<dim>    &dof,
+				      const typename FunctionMap<dim>::type &boundary_functions,
+				      const Quadrature<dim-1>  &q,
+				      std::map<unsigned int,double> &boundary_values)
+{
+//TODO:[?] In VectorTools::project_boundary_values, no condensation of sparsity
+//    structures, matrices and right hand sides or distribution of
+//    solution vectors is performed. This is ok for dim<3 because then
+//    there are no constrained nodes on the boundary, but is not
+//    acceptable for higher dimensions. Fix this.
+
+  Assert (dof.get_fe().n_components() == boundary_functions.begin()->second->n_components,
+	  ExcComponentMismatch());
+  
+  std::vector<unsigned int> dof_to_boundary_mapping;
+  std::set<unsigned char> selected_boundary_components;
+  for (typename FunctionMap<dim>::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);
+  
+				   // set up sparsity structure
+  SparsityPattern sparsity(dof.n_boundary_dofs(boundary_functions),
+			   dof.max_couplings_between_boundary_dofs());
+  DoFTools::make_boundary_sparsity_pattern (dof,
+					    boundary_functions,
+					    dof_to_boundary_mapping,
+					    sparsity);
+
+				   // note: for three or more dimensions, there
+				   // may be constrained nodes on the boundary
+				   // in this case the boundary mass matrix has
+				   // to be condensed and the solution is to
+				   // be distributed afterwards, which is not
+				   // yet implemented. The reason for this is
+				   // that we cannot simply use the @p{condense}
+				   // family of functions, since the matrices
+				   // and vectors do not use the global
+				   // numbering but rather the boundary
+				   // numbering, i.e. the condense function
+				   // needs to use another indirection. There
+				   // should be not many technical problems,
+				   // but it needs to be implemented
+  if (dim<3)
+    sparsity.compress();
+  else
+    Assert (false, ExcNotImplemented());
+
+
+				   // 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, 
+					      mass_matrix, boundary_functions,
+					      rhs, dof_to_boundary_mapping);
+
+				   // same thing as above: if dim>=3 we need
+				   // to consider constraints
+  Assert (dim<3, ExcNotImplemented());
+
+
+  Vector<double> boundary_projection (rhs.size());
+
+  SolverControl           control(1000, 1e-16);
+  PrimitiveVectorMemory<> memory;
+  SolverCG<>              cg(control,memory);
+
+  PreconditionSSOR<> prec;
+  prec.initialize(mass_matrix, 1.2);
+				   // solve
+  cg.solve (mass_matrix, boundary_projection, rhs, prec);
+
+				   // fill in boundary values
+  for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
+    if (dof_to_boundary_mapping[i] != DoFHandler<dim>::invalid_dof_index)
+				       // this dof is on one of the
+				       // interesting boundary parts
+				       //
+				       // remember: @p{i} is the global dof
+				       // number, @p{dof_to_boundary_mapping[i]}
+				       // is the number on the boundary and
+				       // thus in the solution vector
+      boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]);
+}
+
+
+template <int dim>
+void
+VectorTools::project_boundary_values (const DoFHandler<dim>    &dof,
+				      const typename FunctionMap<dim>::type &boundary_functions,
+				      const Quadrature<dim-1>  &q,
+				      std::map<unsigned int,double> &boundary_values)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  project_boundary_values(mapping, dof, boundary_functions, q, boundary_values);
+}
+
+
+
+template <int dim, class InVector, class OutVector>
+void
+VectorTools::integrate_difference (const Mapping<dim>    &mapping,
+				   const DoFHandler<dim> &dof,
+				   const InVector        &fe_function,
+				   const Function<dim>   &exact_solution,
+				   OutVector             &difference,
+				   const Quadrature<dim> &q,
+				   const NormType        &norm,
+				   const Function<dim>   *weight,
+				   const double           exponent_1)
+{
+				   // we mark the "exponent" parameter
+				   // to this function "const" since
+				   // it is strictly incoming, but we
+				   // need to set it to something
+				   // different later on, if
+				   // necessary, so have a read-write
+				   // version of it:
+  double exponent = exponent_1;
+  
+  const unsigned int        n_q_points   = q.n_quadrature_points;
+  const FiniteElement<dim> &fe           = dof.get_fe();
+  const unsigned int        n_components = fe.n_components();
+  const bool                fe_is_system = (n_components != 1);
+
+  if (weight!=0)
+    {
+      Assert ((weight->n_components==1) || (weight->n_components==n_components),
+	      ExcDimensionMismatch(weight->n_components, n_components));
+    }
+
+  difference.reinit (dof.get_tria().n_active_cells());
+  
+  switch (norm)
+    {
+      case L2_norm:
+      case H1_seminorm:
+      case H1_norm:
+	exponent = 2.;
+	break;
+      case L1_norm:
+	exponent = 1.;
+	break;
+      default:
+	break;
+    }
+  
+  UpdateFlags update_flags = UpdateFlags (update_q_points  |
+ 					  update_JxW_values);
+  switch (norm)
+    {
+      case H1_seminorm:
+      case W1p_seminorm:
+      case W1infty_seminorm:
+	update_flags |= UpdateFlags (update_gradients);
+	break;
+      case H1_norm:
+      case W1p_norm:
+      case W1infty_norm:
+	update_flags |= UpdateFlags (update_gradients);
+					 // no break!
+      default:
+	update_flags |= UpdateFlags (update_values);
+	break;
+    }  
+  
+  FEValues<dim> fe_values(mapping, fe, q, update_flags);
+
+  std::vector< Vector<double> >        function_values (n_q_points,
+							Vector<double>(n_components));
+  std::vector<std::vector<Tensor<1,dim> > > function_grads (n_q_points,
+							    std::vector<Tensor<1,dim> >(n_components));
+  std::vector<double> weight_values (n_q_points);
+  std::vector<Vector<double> > weight_vectors (n_q_points, 
+					       Vector<double>(n_components));
+  
+  std::vector<Vector<double> >         psi_values (n_q_points,
+						   Vector<double>(n_components));
+  std::vector<std::vector<Tensor<1,dim> > > psi_grads (n_q_points,
+						       std::vector<Tensor<1,dim> >(n_components));
+  std::vector<double> psi_scalar (n_q_points);
+				   // tmp vector when we use the
+				   // Function<dim> functions for
+				   // scalar functions
+  std::vector<double>         tmp_values (fe_values.n_quadrature_points);
+  std::vector<Tensor<1,dim> > tmp_gradients (fe_values.n_quadrature_points);
+  
+ 				   // loop over all cells
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+						 endc = dof.end();
+  for (unsigned int index=0; cell != endc; ++cell, ++index)
+    {
+      double diff=0;
+ 				       // initialize for this cell
+      fe_values.reinit (cell);
+      
+      if (weight!=0)
+	{
+	  if (weight->n_components>1)
+	    weight->vector_value_list (fe_values.get_quadrature_points(),
+				       weight_vectors);
+	  else
+	    {
+	      weight->value_list (fe_values.get_quadrature_points(),
+				  weight_values);
+	      for (unsigned int k=0;k<n_q_points;++k)
+		weight_vectors[k] = weight_values[k];
+	    }
+	} else {
+	  for (unsigned int k=0;k<n_q_points;++k)
+	    weight_vectors[k] = 1.;
+	}
+      
+      
+      if (update_flags & update_values)
+	{
+					   // first compute the exact solution
+					   // (vectors) at the quadrature points
+					   // try to do this as efficient as
+					   // possible by avoiding a second
+					   // virtual function call in case
+					   // the function really has only
+					   // one component
+	  if (fe_is_system)
+	    exact_solution.vector_value_list (fe_values.get_quadrature_points(),
+					      psi_values);
+	  else
+	    {
+	      exact_solution.value_list (fe_values.get_quadrature_points(),
+					 tmp_values);
+	      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);
+	  for (unsigned int q=0; q<n_q_points; ++q)
+	    psi_values[q] -= function_values[q];
+	}
+
+				       // Do the same for gradients, if required
+      if (update_flags & update_gradients)
+	{
+					   // try to be a little clever
+					   // to avoid recursive virtual
+					   // function calls when calling
+					   // @p{gradient_list} for functions
+					   // that are really scalar
+					   // functions
+	  if (fe_is_system)
+	    exact_solution.vector_gradient_list (fe_values.get_quadrature_points(),
+						 psi_grads);
+	  else
+	    {
+	      exact_solution.gradient_list (fe_values.get_quadrature_points(),
+					    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);
+	  for (unsigned int k=0; k<n_components; ++k)
+	    for (unsigned int q=0; q<n_q_points; ++q)
+	      psi_grads[q][k] -= function_grads[q][k];
+	}
+      
+      switch (norm)
+	{
+	  case mean:
+	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
+					     // Compute values in
+					     // quadrature points
+	    for (unsigned int k=0; k<n_components; ++k)
+	      for (unsigned int q=0; q<n_q_points; ++q)
+		  psi_scalar[q] += psi_values[q](k)
+				   * weight_vectors[q](k);
+
+					     // Integrate
+	    diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
+				       fe_values.get_JxW_values().begin(),
+				       0.0);
+	    break;
+	  case Lp_norm:
+	  case L1_norm:
+	  case W1p_norm:
+	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
+					     // Compute values in
+					     // quadrature points
+	    for (unsigned int k=0; k<n_components; ++k)
+	      for (unsigned int q=0; q<n_q_points; ++q)
+		  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(),
+				       0.0);
+					     // Compute the root only,
+					     // if no derivative
+					     // values are added later
+	    if (!(update_flags & update_gradients))
+	      diff = std::pow(diff, 1./exponent);
+	    break;
+	  case L2_norm:
+	  case H1_norm:
+	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
+					     // Compute values in
+					     // quadrature points
+	    for (unsigned int k=0; k<n_components; ++k)
+	      for (unsigned int q=0; q<n_q_points; ++q)
+		  psi_scalar[q] += psi_values[q](k)*psi_values[q](k)
+				   * weight_vectors[q](k);
+
+					     // Integrate
+	    diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
+				       fe_values.get_JxW_values().begin(),
+				       0.0);
+					     // Compute the root only,
+					     // if no derivative
+					     // values are added later
+	    if (norm == L2_norm)
+	      diff=std::sqrt(diff);
+	    break;
+	  case Linfty_norm:
+	  case W1infty_norm:
+	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
+	    for (unsigned int k=0; k<n_components; ++k)
+	      for (unsigned int q=0; q<n_q_points; ++q)
+		{
+		  double newval = std::fabs(psi_values[q](k))
+				  * weight_vectors[q](k);
+		  if (psi_scalar[q]<newval)
+		    psi_scalar[q] = newval;
+		}
+					     // Maximum on one cell
+	    diff = *std::max_element (psi_scalar.begin(), psi_scalar.end());
+	    break;
+	  case H1_seminorm:
+	  case W1p_seminorm:
+	    break;
+	  default:
+	    Assert (false, ExcNotImplemented());
+	    break;
+	}
+
+      switch (norm)
+	{
+	  case W1p_seminorm:
+	  case W1p_norm:
+	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
+	    for (unsigned int k=0; k<n_components; ++k)
+	      for (unsigned int q=0; q<n_q_points; ++q)
+		psi_scalar[q] += std::pow(sqr_point(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);
+	    diff = std::pow(diff, 1./exponent);
+	    break;
+	  case H1_seminorm:
+	  case H1_norm:
+					     // take square of integrand
+	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
+	    for (unsigned int k=0; k<n_components; ++k)
+	      for (unsigned int q=0; q<n_q_points; ++q)
+		psi_scalar[q] += sqr_point(psi_grads[q][k])
+				 * weight_vectors[q](k);
+
+ 					     // add seminorm to L_2 norm or
+ 					     // to zero
+	    diff += std::inner_product (psi_scalar.begin(), psi_scalar.end(),
+					fe_values.get_JxW_values().begin(),
+					0.0);
+	    diff = std::sqrt(diff);
+	    break;
+	  case W1infty_seminorm:
+	  case W1infty_norm:
+	    Assert(false, ExcNotImplemented());
+	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
+	    for (unsigned int k=0; k<n_components; ++k)
+	      for (unsigned int q=0; q<n_q_points; ++q)
+		{
+		  double t = 0.;
+		  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);
+		}
+
+	    for (unsigned int i=0;i<psi_scalar.size();++i)
+	      diff = std::max (diff, psi_scalar[i]);
+	    break;
+	  default:
+	    break;
+ 	}
+				       // append result of this cell
+ 				       // to the end of the vector
+      difference(index) = diff;
+    }
+}
+
+
+template <int dim, class InVector, class OutVector>
+void
+VectorTools::integrate_difference (const DoFHandler<dim>    &dof,
+				   const InVector           &fe_function,
+				   const Function<dim>      &exact_solution,
+				   OutVector                &difference,
+				   const Quadrature<dim>    &q,
+				   const NormType           &norm,
+				   const Function<dim>      *weight,
+				   const double              exponent)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  integrate_difference(mapping, dof, fe_function, exact_solution,
+		       difference, q, norm, weight, exponent);
+}
+
+
+
+template <int dim, class InVector>
+void
+VectorTools::point_difference (const DoFHandler<dim> &dof,
+			       const InVector        &fe_function,
+			       const Function<dim>   &exact_function,
+			       Vector<double>        &difference,
+			       const Point<dim>      &point)
+{
+  static const MappingQ1<dim> mapping;
+  const FiniteElement<dim>& fe = dof.get_fe();
+
+  Assert(difference.size() == fe.n_components(),
+	 ExcDimensionMismatch(difference.size(), fe.n_components()));
+
+                                   // first find the cell in which this point
+                                   // is, initialize a quadrature rule with
+                                   // it, and then a FEValues object
+  const typename DoFHandler<dim>::active_cell_iterator
+    cell = GridTools::find_active_cell_around_point (dof, point);
+
+  const Point<dim> unit_point
+    = mapping.transform_real_to_unit_cell(cell, point);
+  Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
+          ExcInternalError());
+
+  const Quadrature<dim> quadrature (std::vector<Point<dim> > (1, unit_point),
+                                    std::vector<double> (1, 1.));
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+  fe_values.reinit(cell);
+
+                                   // then use this to get at the values of
+                                   // the given fe_function at this point
+  std::vector<Vector<double> > u_value(1, Vector<double> (fe.n_components()));
+  fe_values.get_function_values(fe_function, u_value);
+
+  if (fe.n_components() == 1)
+    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);
+}
+
+
+
+template <int dim, class InVector>
+double
+VectorTools::compute_mean_value (const Mapping<dim>    &mapping,
+				 const DoFHandler<dim> &dof,
+				 const Quadrature<dim> &quadrature,
+				 const InVector        &v,
+				 const unsigned int     component)
+{
+  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));
+
+  typename DoFHandler<dim>::active_cell_iterator c;
+  std::vector<Vector<double> > values(quadrature.n_quadrature_points,
+				      Vector<double> (dof.get_fe().n_components()));
+  
+  double mean = 0.;
+  double area = 0.;
+				   // Compute mean value
+  for (c = dof.begin_active(); c != dof.end(); ++c)
+    {
+      fe.reinit (c);
+      fe.get_function_values(v, values);
+      for (unsigned int k=0; k< quadrature.n_quadrature_points; ++k)
+	{
+	  mean += fe.JxW(k) * values[k](component);
+	  area += fe.JxW(k);
+	};
+    };
+  
+  return (mean/area);
+}
+
+
+template <int dim, class InVector>
+double
+VectorTools::compute_mean_value (const DoFHandler<dim> &dof,
+				 const Quadrature<dim> &quadrature,
+				 const InVector        &v,
+				 const unsigned int     component)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  return compute_mean_value(mapping, dof, quadrature, v, component);
+}
+
+#endif
diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc
index 2bbad617ac..3dbc4d1779 100644
--- a/deal.II/deal.II/source/numerics/vectors.cc
+++ b/deal.II/deal.II/source/numerics/vectors.cc
@@ -11,1705 +11,25 @@
 //
 //---------------------------------------------------------------------------
 
+#include<numerics/vectors.templates.h>
 
-#include <base/function.h>
-#include <base/quadrature.h>
-#include <lac/vector.h>
-#include <lac/block_vector.h>
-#include <lac/sparse_matrix.h>
-#include <lac/precondition.h>
-#include <lac/solver_cg.h>
-#include <lac/vector_memory.h>
-#include <grid/tria_iterator.h>
-#include <grid/grid_tools.h>
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-#include <fe/fe.h>
-#include <fe/fe_values.h>
-#include <fe/mapping_q1.h>
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-
-#include <numeric>
-#include <algorithm>
-#include <vector>
-#include <cmath>
-
-//TODO:[GK] Move templates containing vector arguments to vectors.templates.h
-
-
-
-template <int dim>
-inline double sqr_point (const Tensor<1,dim> &p)
-{
-  return p * p;
-}
-
-
-
-
-template <int dim, class VECTOR>
-void VectorTools::interpolate (const Mapping<dim>    &mapping,
-			       const DoFHandler<dim> &dof,
-			       const Function<dim>   &function,
-			       VECTOR                &vec)
-{
-  Assert (dof.get_fe().n_components() == function.n_components,
-	  ExcComponentMismatch());
-  
-  const FiniteElement<dim> &fe           = dof.get_fe();
-  const unsigned int        n_components = fe.n_components();
-  const bool                fe_is_system = (n_components != 1);
-  
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
-						 endc = dof.end();
-
-				   // For FESystems many of the
-				   // unit_support_points will
-				   // appear multiply, as a point
-				   // may be unit_support_point
-				   // for several of the components
-				   // of the system.
-				   // The following is rather
-				   // complicated as it is
-				   // avoided to evaluate
-				   // the vectorfunction multiply at
-				   // the same point on a cell.
-  const std::vector<Point<dim> > &
-    unit_support_points = fe.get_unit_support_points();
-  Assert (unit_support_points.size() != 0,
-	  ExcNonInterpolatingFE());
-
-				   // Find the support points 
-				   // on a cell that
-				   // are multiply mentioned in 
-				   // @p{unit_support_points}.
-				   // Mark the first representative
-				   // of each multiply mentioned
-				   // support point by appending its
-				   // dof index to @p{dofs_of_rep_points}.
-				   // Each multiple point gets to know
-				   // the dof index of its representative
-				   // point by the @p{dof_to_rep_dof_table}.
-
-				   // the following vector collects all dofs i,
-				   // 0<=i<fe.dofs_per_cell, for that
-				   // 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
-				   // is called rep index.
-  std::vector<unsigned int> dofs_of_rep_points;
-				   // the following table converts a dof i
-				   // to the rep index.
-  std::vector<unsigned int> dof_to_rep_index_table;
-  unsigned int n_rep_points=0;
-  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-    {
-      bool representative=true;
-				       // the following loop is looped
-				       // the other way round to get
-				       // the minimal effort of
-				       // O(fe.dofs_per_cell) for multiple
-				       // support points that are placed
-				       // one after the other.
-      for (unsigned int j=dofs_of_rep_points.size(); j>0; --j)
-	if (unit_support_points[i] 
-	    == unit_support_points[dofs_of_rep_points[j-1]])
-	  {
-	    dof_to_rep_index_table.push_back(j-1);
-	    representative=false;
-	    break;
-	  }
-      
-      if (representative)
-	{
-					   // rep_index=dofs_of_rep_points.size()
-	  dof_to_rep_index_table.push_back(dofs_of_rep_points.size());
-					   // dofs_of_rep_points[rep_index]=i
-	  dofs_of_rep_points.push_back(i);
-	  ++n_rep_points;
-	}
-    }
-  Assert(dofs_of_rep_points.size()==n_rep_points, ExcInternalError());
-  Assert(dof_to_rep_index_table.size()==fe.dofs_per_cell, ExcInternalError());
-
-  std::vector<unsigned int> dofs_on_cell (fe.dofs_per_cell);
-  std::vector<Point<dim> >  rep_points (n_rep_points);
-
-				   // get space for the values of the
-				   // function at the rep support points.
-				   //
-				   // have two versions, one for system fe
-				   // and one for scalar ones, to take the
-				   // more efficient one respectively
-  std::vector<double>          function_values_scalar (n_rep_points);
-  std::vector<Vector<double> > function_values_system (n_rep_points,
-						  Vector<double>(fe.n_components()));
-
-				   // Make a quadrature rule from support points
-				   // to feed it into FEValues
-  Quadrature<dim> support_quadrature(unit_support_points);
-
-				   // Transformed support points are computed by
-				   // FEValues
-  FEValues<dim> fe_values (mapping, fe, support_quadrature, update_q_points);
-  
-  for (; cell!=endc; ++cell)
-    {
-				       // for each cell:
-				       // get location of finite element
-				       // support_points
-      fe_values.reinit(cell);
-      const std::vector<Point<dim> >& support_points =
-	fe_values.get_quadrature_points();
-      
-				       // pick out the representative
-				       // support points
-      for (unsigned int j=0; j<dofs_of_rep_points.size(); ++j)
-	rep_points[j]=support_points[dofs_of_rep_points[j]];
-
-				       // get indices of the dofs on this cell
-      cell->get_dof_indices (dofs_on_cell);
-
-
-      if (fe_is_system)
-	{
-					   // get function values at
-					   // these points. Here: get
-					   // all components
-	  function.vector_value_list (rep_points, function_values_system);
-					   // distribute the function
-					   // values to the global
-					   // vector
-	  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-	    {
-	      const unsigned int component
-		= fe.system_to_component_index(i).first;
-	      const unsigned int rep_dof=dof_to_rep_index_table[i];
-	      vec(dofs_on_cell[i])
-		= function_values_system[rep_dof](component);
-	    };
-	}
-      
-      else
-	{
-					   // get first component only,
-					   // which is the only component
-					   // in the function anyway
-	  function.value_list (rep_points, function_values_scalar, 0);
-					   // distribute the function
-					   // values to the global
-					   // vector
-	  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-	    vec(dofs_on_cell[i]) 
-	      = function_values_scalar[dof_to_rep_index_table[i]];
-	};
-    }
-}
-
-
-template <int dim, class VECTOR>
-void VectorTools::interpolate (const DoFHandler<dim> &dof,
-			       const Function<dim>   &function,
-			       VECTOR                &vec)
-{
-  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  static const MappingQ1<dim> mapping;
-  interpolate(mapping, dof, function, vec);
-}
-
-
-
-
-template <int dim, class InVector, class OutVector>
-void
-VectorTools::interpolate (const DoFHandler<dim>           &dof_1,
-			  const DoFHandler<dim>           &dof_2,
-			  const FullMatrix<double>        &transfer,
-			  const InVector                  &data_1,
-			  OutVector                       &data_2)
-{
-  Vector<double> cell_data_1(dof_1.get_fe().dofs_per_cell);
-  Vector<double> cell_data_2(dof_2.get_fe().dofs_per_cell);
-
-  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>::active_cell_iterator h = dof_1.begin_active();
-  typename DoFHandler<dim>::active_cell_iterator l = dof_2.begin_active();
-  const typename DoFHandler<dim>::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) 
-      {
-	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());	
-	++touch_count[local_dof_indices[j]];
-      };
-  };
-
-				   // compute the mean value of the
-				   // sum which we have placed in each
-				   // entry of the output vector
-  for (unsigned int i=0; i<dof_2.n_dofs(); ++i)
-    {
-      Assert (touch_count[i] != 0,
-	      ExcInternalError());
-      
-      data_2(i) /= touch_count[i];
-    };
-}
-
-
-#if deal_II_dimension == 1
-
-void VectorTools::project (const Mapping<1>       &,
-			   const DoFHandler<1>    &,
-			   const ConstraintMatrix &,
-			   const Quadrature<1>    &,
-			   const Function<1>      &,
-			   Vector<double>         &,
-			   const bool              ,
-			   const Quadrature<0>    &,
-			   const bool              )
-{
-				   // this function should easily be implemented
-				   // using the template below. However some
-				   // changes have to be made since faces don't
-				   // exist in 1D. Maybe integrate the creation of
-				   // zero boundary values into the
-				   // project_boundary_values function?
-  Assert (false, ExcNotImplemented());
-}
-
-
-#endif
-
-
-template <int dim>
-void VectorTools::project (const Mapping<dim>       &mapping,
-			   const DoFHandler<dim>    &dof,
-			   const ConstraintMatrix   &constraints,
-			   const Quadrature<dim>    &quadrature,
-			   const Function<dim>      &function,
-			   Vector<double>           &vec,
-			   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,
-	  ExcInvalidFE());
-  
-  const FiniteElement<dim> &fe = dof.get_fe();
-
-				   // make up boundary values
-  std::map<unsigned int,double> boundary_values;
-
-  if (enforce_zero_boundary == true) 
-				     // no need to project boundary
-				     // values, but enforce
-				     // homogeneous boundary values
-				     // anyway
-    {
-				       // loop over all boundary faces
-				       // to get all dof indices of
-				       // dofs on the boundary. note
-				       // that in 3d there are cases
-				       // where a face is not at the
-				       // boundary, yet one of its
-				       // lines is, and we should
-				       // consider the degrees of
-				       // freedom on it as boundary
-				       // nodes. likewise, in 2d and
-				       // 3d there are cases where a
-				       // cell is only at the boundary
-				       // by one vertex. nevertheless,
-				       // since we do not support
-				       // boundaries with dimension
-				       // less or equal to dim-2, each
-				       // such boundary dof is also
-				       // found from some other face
-				       // that is actually wholly on
-				       // the boundary, not only by
-				       // one line or one vertex
-      typename DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
-						     endf = dof.end_face();
-      std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
-      for (; face!=endf; ++face)
-	if (face->at_boundary())
-	  {
-	    face->get_dof_indices (face_dof_indices);
-	    for (unsigned int i=0; i<fe.dofs_per_face; ++i)
-					       // enter zero boundary values
-					       // for all boundary nodes
-					       //
-					       // we need not care about
-					       // vector valued elements here,
-					       // since we set all components
-	      boundary_values[face_dof_indices[i]] = 0.;
-	  };
-    }
-  else
-				     // no homogeneous boundary values
-    if (project_to_boundary_first == true)
-				       // boundary projection required
-      {
-					 // set up a list of boundary functions for
-					 // the different boundary parts. We want the
-					 // @p{function} to hold on all parts of the
-					 // boundary
-	typename FunctionMap<dim>::type boundary_functions;
-	for (unsigned char c=0; c<255; ++c)
-	  boundary_functions[c] = &function;
-	project_boundary_values (dof, boundary_functions, q_boundary,
-				 boundary_values);
-      };
-
-
-				   // set up mass matrix and right hand side
-  vec.reinit (dof.n_dofs());
-  SparsityPattern sparsity(dof.n_dofs(),
-			   dof.n_dofs(),
-			   dof.max_couplings_between_dofs());
-  DoFTools::make_sparsity_pattern (dof, sparsity);
-  constraints.condense (sparsity);
-  
-  SparseMatrix<double> mass_matrix (sparsity);
-  Vector<double> tmp (mass_matrix.n());
-
-  MatrixCreator::create_mass_matrix (mapping, dof, quadrature, mass_matrix);
-  
-  VectorTools::create_right_hand_side (mapping, dof, quadrature, function, tmp);
-
-  constraints.condense (mass_matrix);
-  constraints.condense (tmp);
-  if (boundary_values.size() != 0)
-    MatrixTools::apply_boundary_values (boundary_values,
-					mass_matrix, vec, tmp,
-					true);
-
-  SolverControl           control(1000,1e-16);
-  PrimitiveVectorMemory<> memory;
-  SolverCG<>              cg(control,memory);
-
-  PreconditionSSOR<> prec;
-  prec.initialize(mass_matrix, 1.2);
-				   // solve
-  cg.solve (mass_matrix, vec, tmp, prec);
-  
-				   // distribute solution
-  constraints.distribute (vec);
-}
-
-
-template <int dim>
-void VectorTools::project (const DoFHandler<dim>    &dof,
-			   const ConstraintMatrix   &constraints,
-			   const Quadrature<dim>    &quadrature,
-			   const Function<dim>      &function,
-			   Vector<double>           &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"));
-  static const MappingQ1<dim> mapping;
-  project(mapping, dof, constraints, quadrature, function, vec,
-	  enforce_zero_boundary, q_boundary, project_to_boundary_first);
-}
-
-
-
-
-template <int dim>
-void VectorTools::create_right_hand_side (const Mapping<dim>    &mapping,
-					  const DoFHandler<dim> &dof_handler,
-					  const Quadrature<dim> &quadrature,
-					  const Function<dim>   &rhs_function,
-					  Vector<double>        &rhs_vector)
-{
-  const FiniteElement<dim> &fe  = dof_handler.get_fe();
-  Assert (fe.n_components() == rhs_function.n_components,
-	  ExcComponentMismatch());
-  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_q_points |
-					 update_JxW_values);
-  FEValues<dim> fe_values (mapping, fe, quadrature, update_flags);
-
-  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>::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) 
-	{
-	  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) 
-	      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));
-      
-      // Use the faster code if the FiniteElement is primitive
-      if (fe.is_primitive ())
-	{
-	  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;
-	      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];
-		  }
-	      
-	      cell->get_dof_indices (dofs);
-	      
-	      for (unsigned int i=0; i<dofs_per_cell; ++i)
-		rhs_vector(dofs[i]) += cell_vector(i);
-	    }
-	}
-      else
-	// Otherwise do it the way proposed for vector valued elements
-	{
-	  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;
-	      for (unsigned int point=0; point<n_q_points; ++point)
-		for (unsigned int i=0; i<dofs_per_cell; ++i)
-		  for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-		    if (fe.get_nonzero_components(i)[comp_i])
-		      {
-			cell_vector(i) += rhs_values[point](comp_i) *
-			                  fe_values.shape_value_component(i,point,comp_i) *
-			                  weights[point];
-		      }
-	      
-	      cell->get_dof_indices (dofs);
-	      
-	      for (unsigned int i=0; i<dofs_per_cell; ++i)
-		rhs_vector(dofs[i]) += cell_vector(i);
-	    }
-	}
-    }
-}
-
-
-template <int dim>
-void VectorTools::create_right_hand_side (const DoFHandler<dim>    &dof_handler,
-					  const Quadrature<dim>    &quadrature,
-					  const Function<dim>      &rhs_function,
-					  Vector<double>           &rhs_vector)
-{
-  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  static const MappingQ1<dim> mapping;
-  create_right_hand_side(mapping, dof_handler, quadrature,
-			 rhs_function, rhs_vector);
-}
-
-
-
-#if deal_II_dimension != 1
-
-template <int dim>
-void
-VectorTools::create_boundary_right_hand_side (const Mapping<dim>      &mapping,
-					      const DoFHandler<dim>   &dof_handler,
-					      const Quadrature<dim-1> &quadrature,
-					      const Function<dim>     &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,
-	  ExcComponentMismatch());
-  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_q_points |
-					 update_JxW_values);
-  FEFaceValues<dim> fe_values (mapping, fe, quadrature, update_flags);
-
-  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>::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 () &&
-	      (boundary_indicators.find (cell->face(face)->boundary_indicator())
-	       !=
-	       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) 
-		  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 (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())
-	       !=
-	       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];
-		      }
-		}
-	      else
-		{
-		  // And the full featured code, if vector valued FEs are used
-		  for (unsigned int point=0; point<n_q_points; ++point)
-		    for (unsigned int i=0; i<dofs_per_cell; ++i)
-		      for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-			if (fe.get_nonzero_components(i)[comp_i])
-			  {
-			    cell_vector(i) += rhs_values[point](comp_i) *
-			                      fe_values.shape_value_component(i,point,comp_i) *
-			                      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
-
-void
-VectorTools::create_boundary_right_hand_side (const Mapping<1>    &,
-					      const DoFHandler<1> &,
-					      const Quadrature<0> &,
-					      const Function<1>   &,
-					      Vector<double>      &,
-					      const std::set<unsigned char> &)
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-#endif
-
-template <int dim>
-void
-VectorTools::create_boundary_right_hand_side (const DoFHandler<dim>   &dof_handler,
-					      const Quadrature<dim-1> &quadrature,
-					      const Function<dim>     &rhs_function,
-					      Vector<double>          &rhs_vector,
-					      const std::set<unsigned char> &boundary_indicators)
-{
-  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  static const MappingQ1<dim> mapping;
-  create_boundary_right_hand_side(mapping, dof_handler, quadrature,
-				  rhs_function, rhs_vector,
-				  boundary_indicators);
-}
-
-
-
-#if deal_II_dimension == 1
-
-void
-VectorTools::interpolate_boundary_values (const Mapping<1>         &,
-					  const DoFHandler<1>      &dof,
-					  const unsigned char       boundary_component,
-					  const Function<1>        &boundary_function,
-					  std::map<unsigned int,double> &boundary_values,
-					  const std::vector<bool>       &component_mask_)
-{
-  Assert (boundary_component != 255,
-	  ExcInvalidBoundaryIndicator());
-
-  const FiniteElement<1> &fe = dof.get_fe();
-  Assert (fe.n_components() == boundary_function.n_components,
-	  ExcComponentMismatch());
-
-				   // set the component mask to either
-				   // the original value or a vector
-				   // of @p{true}s
-  const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
-					  std::vector<bool> (fe.n_components(), true) :
-					  component_mask_);
-  Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
-	  ExcComponentMismatch());
-  
-				   // check whether boundary values at
-				   // the left or right boundary of
-				   // the line are
-				   // requested. @p{direction} denotes
-				   // the neighboring direction in
-				   // which we seek the boundary,
-				   // i.e. 0 is left boundary and 1 is
-				   // 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
-				   // to the children
-  DoFHandler<1>::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);
-
-				   // now set the value of the
-				   // outermost degree of
-				   // freedom. setting also
-				   // creates the entry in the map
-				   // if it did not exist
-				   // beforehand
-				   //
-				   // save some time by requesting
-				   // values only once for each point,
-				   // irrespective of the number of
-				   // components of the function
-  Vector<double> function_values (fe.n_components());
-  if (fe.n_components() == 1)
-    function_values(0)
-      = boundary_function.value (outermost_cell->vertex(direction));
-  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)]
-	= function_values(fe.face_system_to_component_index(i).first);
-}
-
-
-
-void
-VectorTools::interpolate_boundary_values (const Mapping<1>              &mapping,
-					  const DoFHandler<1>           &dof,
-					  const FunctionMap<1>::type    &function_map,
-					  std::map<unsigned int,double> &boundary_values,
-					  const std::vector<bool>       &component_mask)
-{
-  for (FunctionMap<1>::type::const_iterator i=function_map.begin();
-       i!=function_map.end(); ++i)
-    interpolate_boundary_values (mapping, dof, i->first, *i->second,
-				 boundary_values, component_mask);
-}
-
-
-#endif
-
-
-template <int dim>
-void
-VectorTools::
-interpolate_boundary_values (const Mapping<dim>            &mapping,
-                             const DoFHandler<dim>         &dof,
-                             const typename FunctionMap<dim>::type &function_map,
-                             std::map<unsigned int,double> &boundary_values,
-                             const std::vector<bool>       &component_mask_)
-{
-				   // if for whatever reason we were
-				   // passed an empty map, return
-				   // immediately
-  if (function_map.size() == 0)
-    return;
-  
-  Assert (function_map.find(255) == function_map.end(),
-	  ExcInvalidBoundaryIndicator());
-
-  const FiniteElement<dim> &fe           = dof.get_fe();
-  const unsigned int        n_components = fe.n_components();
-  const bool                fe_is_system = (n_components != 1);
-
-  for (typename FunctionMap<dim>::type::const_iterator i=function_map.begin();
-       i!=function_map.end(); ++i)
-    Assert (n_components == i->second->n_components,
-	    ExcInvalidFE());
-
-				   // set the component mask to either
-				   // the original value or a vector
-				   // of @p{true}s
-  const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
-					  std::vector<bool> (fe.n_components(), true) :
-					  component_mask_);
-  Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
-	  ExcComponentMismatch());
-
-				   // field to store the indices
-  std::vector<unsigned int> face_dofs (fe.dofs_per_face,
-				       DoFHandler<dim>::invalid_dof_index);
-  std::vector<Point<dim> >  dof_locations (face_dofs.size(), Point<dim>());
-  
-				   // array to store the values of
-				   // the boundary function at the
-				   // boundary points. have to arrays
-				   // for scalar and vector functions
-				   // to use the more efficient one
-				   // respectively
-  std::vector<double>          dof_values_scalar (fe.dofs_per_face);
-  std::vector<Vector<double> > dof_values_system (fe.dofs_per_face,
-						  Vector<double>(fe.n_components()));
-
-				   // next generate a quadrature rule
-				   // on the face from the unit
-				   // support points. this wil be used
-				   // to obtain the quadrature points
-				   // on the real cell's face
-  std::vector<Point<dim-1> >
-    unit_support_points = fe.get_unit_face_support_points();
-  
-				   // check whether there are support
-				   // points on the face. if not, then
-				   // we should try a more clever
-				   // way. the idea is that a finite
-				   // element may not offer support
-				   // points for all its shape
-				   // functions, but maybe only
-				   // some. if it offers support
-				   // points for the components we are
-				   // interested in in this function,
-				   // then that's fine. if not, the
-				   // function we call in the finite
-				   // element will raise an
-				   // exception. the support points
-				   // for the other shape functions
-				   // are left uninitialized (well,
-				   // initialized by the default
-				   // constructor), since we don't
-				   // need them anyway.
-  if (unit_support_points.size() == 0)
-    {
-      unit_support_points.resize (fe.dofs_per_face);
-      for (unsigned int i=0; i<fe.dofs_per_face; ++i)
-        if (component_mask[fe.face_system_to_component_index(i).first]
-            == true)
-          unit_support_points[i] = fe.unit_face_support_point(i);
-    };
-
-  Quadrature<dim-1> aux_quad (unit_support_points);
-  FEFaceValues<dim> fe_values (mapping, fe, aux_quad, update_q_points);
-
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
-						 endc = dof.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
-	 ++face_no)
-      {
-	typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
-	const unsigned char boundary_component = face->boundary_indicator();
-	if (function_map.find(boundary_component) != function_map.end()) 
-					   // face is of the right component
-	  {
-					     // get indices, physical location and
-					     // boundary values of dofs on this
-					     // face
-	    face->get_dof_indices (face_dofs);
-	    fe_values.reinit(cell, face_no);
-	    const std::vector<Point<dim> > &dof_locations = fe_values.get_quadrature_points ();
-	    
-	    if (fe_is_system)
-	      {
-		function_map.find(boundary_component)->second
-                  ->vector_value_list (dof_locations, dof_values_system);
-		
-						 // enter those dofs
-						 // into the list that
-						 // match the
-						 // component
-						 // signature. avoid
-						 // the usual
-						 // complication that
-						 // we can't just use
-						 // *_system_to_component_index
-						 // for non-primitive
-						 // FEs
-		for (unsigned int i=0; i<face_dofs.size(); ++i)
-                  {
-                    unsigned int component;
-                    if (fe.is_primitive())
-                      component = fe.face_system_to_component_index(i).first;
-                    else
-                      {
-                                                         // non-primitive
-                                                         // case. make
-                                                         // sure that
-                                                         // this
-                                                         // particular
-                                                         // shape
-                                                         // function
-                                                         // _is_
-                                                         // primitive,
-                                                         // and get at
-                                                         // it's
-                                                         // component. use
-                                                         // usual
-                                                         // trick to
-                                                         // transfer
-                                                         // face dof
-                                                         // index to
-                                                         // cell dof
-                                                         // index
-                        const unsigned int cell_i
-                          = (dim == 1 ?
-                             i
-                             :
-                             (dim == 2 ?
-                              (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
-                              :
-                              (dim == 3 ?
-                               (i<4*fe.dofs_per_vertex ?
-                                i
-                                :
-                                (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
-                                 i+4*fe.dofs_per_vertex
-                                 :
-                                 i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
-                               :
-                               deal_II_numbers::invalid_unsigned_int)));
-                        Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
-
-                                                         // make sure
-                                                         // that if
-                                                         // this is
-                                                         // not a
-                                                         // primitive
-                                                         // shape function,
-                                                         // then all
-                                                         // the
-                                                         // corresponding
-                                                         // components
-                                                         // in the
-                                                         // mask are
-                                                         // not set
-                        if (!fe.is_primitive(cell_i))
-                          for (unsigned int c=0; c<n_components; ++c)
-                            if (fe.get_nonzero_components(cell_i)[c])
-                              Assert (component_mask[c] == false,
-                                      ExcFENotPrimitive());
-
-                                                         // let's pick
-                                                         // the first
-                                                         // of
-                                                         // possibly
-                                                         // more than
-                                                         // one
-                                                         // non-zero
-                                                         // components. if
-                                                         // shape
-                                                         // function
-                                                         // is
-                                                         // non-primitive,
-                                                         // then we
-                                                         // will
-                                                         // ignore the
-                                                         // result in
-                                                         // the
-                                                         // following
-                                                         // anyway,
-                                                         // otherwise
-                                                         // there's
-                                                         // only one
-                                                         // non-zero
-                                                         // component
-                                                         // which we
-                                                         // will use
-                        component = (std::find (fe.get_nonzero_components(cell_i).begin(),
-                                                fe.get_nonzero_components(cell_i).end(),
-                                                true)
-                                     -
-                                     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,
-					       // so save some computations
-	      {
-						 // get only the one component that
-						 // this function has
-		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];
-	      }
-	  }
-      }
-}
-
-
-
-template <int dim>
-void
-VectorTools::interpolate_boundary_values (const Mapping<dim>            &mapping,
-					  const DoFHandler<dim>         &dof,
-					  const unsigned char            boundary_component,
-					  const Function<dim>           &boundary_function,
-					  std::map<unsigned int,double> &boundary_values,
-					  const std::vector<bool>       &component_mask)
-{
-  typename FunctionMap<dim>::type function_map;
-  function_map[boundary_component] = &boundary_function;
-  interpolate_boundary_values (mapping, dof, function_map, boundary_values,
-			       component_mask);
-}
-
-
-  
-template <int dim>
-void
-VectorTools::interpolate_boundary_values (const DoFHandler<dim>         &dof,
-					  const unsigned char            boundary_component,
-					  const Function<dim>           &boundary_function,
-					  std::map<unsigned int,double> &boundary_values,
-					  const std::vector<bool>       &component_mask)
-{
-  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  static const MappingQ1<dim> mapping;
-  interpolate_boundary_values(mapping, dof, boundary_component,
-			      boundary_function, boundary_values, component_mask);
-}
-
-
-
-template <int dim>
-void
-VectorTools::interpolate_boundary_values (const DoFHandler<dim>         &dof,
-					  const typename FunctionMap<dim>::type &function_map,
-					  std::map<unsigned int,double> &boundary_values,
-					  const std::vector<bool>       &component_mask)
-{
-  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  static const MappingQ1<dim> mapping;
-  interpolate_boundary_values(mapping, dof, function_map,
-			      boundary_values, component_mask);
-}
-
-
-#if deal_II_dimension == 1
-
-void
-VectorTools::project_boundary_values (const Mapping<1>       &mapping,
-				      const DoFHandler<1>    &dof,
-				      const FunctionMap<1>::type &boundary_functions,
-				      const Quadrature<0>  &,
-				      std::map<unsigned int,double> &boundary_values)
-{
-				   // projection in 1d is equivalent
-				   // to interpolation
-  interpolate_boundary_values (mapping, dof, boundary_functions,
-			       boundary_values, std::vector<bool>());
-}
-
-#endif
-
-
-template <int dim>
-void
-VectorTools::project_boundary_values (const Mapping<dim>       &mapping,
-				      const DoFHandler<dim>    &dof,
-				      const typename FunctionMap<dim>::type &boundary_functions,
-				      const Quadrature<dim-1>  &q,
-				      std::map<unsigned int,double> &boundary_values)
-{
-//TODO:[?] In VectorTools::project_boundary_values, no condensation of sparsity
-//    structures, matrices and right hand sides or distribution of
-//    solution vectors is performed. This is ok for dim<3 because then
-//    there are no constrained nodes on the boundary, but is not
-//    acceptable for higher dimensions. Fix this.
-
-  Assert (dof.get_fe().n_components() == boundary_functions.begin()->second->n_components,
-	  ExcComponentMismatch());
-  
-  std::vector<unsigned int> dof_to_boundary_mapping;
-  std::set<unsigned char> selected_boundary_components;
-  for (typename FunctionMap<dim>::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);
-  
-				   // set up sparsity structure
-  SparsityPattern sparsity(dof.n_boundary_dofs(boundary_functions),
-			   dof.max_couplings_between_boundary_dofs());
-  DoFTools::make_boundary_sparsity_pattern (dof,
-					    boundary_functions,
-					    dof_to_boundary_mapping,
-					    sparsity);
-
-				   // note: for three or more dimensions, there
-				   // may be constrained nodes on the boundary
-				   // in this case the boundary mass matrix has
-				   // to be condensed and the solution is to
-				   // be distributed afterwards, which is not
-				   // yet implemented. The reason for this is
-				   // that we cannot simply use the @p{condense}
-				   // family of functions, since the matrices
-				   // and vectors do not use the global
-				   // numbering but rather the boundary
-				   // numbering, i.e. the condense function
-				   // needs to use another indirection. There
-				   // should be not many technical problems,
-				   // but it needs to be implemented
-  if (dim<3)
-    sparsity.compress();
-  else
-    Assert (false, ExcNotImplemented());
-
-
-				   // 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, 
-					      mass_matrix, boundary_functions,
-					      rhs, dof_to_boundary_mapping);
-
-				   // same thing as above: if dim>=3 we need
-				   // to consider constraints
-  Assert (dim<3, ExcNotImplemented());
-
-
-  Vector<double> boundary_projection (rhs.size());
-
-  SolverControl           control(1000, 1e-16);
-  PrimitiveVectorMemory<> memory;
-  SolverCG<>              cg(control,memory);
-
-  PreconditionSSOR<> prec;
-  prec.initialize(mass_matrix, 1.2);
-				   // solve
-  cg.solve (mass_matrix, boundary_projection, rhs, prec);
-
-				   // fill in boundary values
-  for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
-    if (dof_to_boundary_mapping[i] != DoFHandler<dim>::invalid_dof_index)
-				       // this dof is on one of the
-				       // interesting boundary parts
-				       //
-				       // remember: @p{i} is the global dof
-				       // number, @p{dof_to_boundary_mapping[i]}
-				       // is the number on the boundary and
-				       // thus in the solution vector
-      boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]);
-}
-
-
-template <int dim>
-void
-VectorTools::project_boundary_values (const DoFHandler<dim>    &dof,
-				      const typename FunctionMap<dim>::type &boundary_functions,
-				      const Quadrature<dim-1>  &q,
-				      std::map<unsigned int,double> &boundary_values)
-{
-  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  static const MappingQ1<dim> mapping;
-  project_boundary_values(mapping, dof, boundary_functions, q, boundary_values);
-}
-
-
-
-template <int dim, class InVector, class OutVector>
-void
-VectorTools::integrate_difference (const Mapping<dim>    &mapping,
-				   const DoFHandler<dim> &dof,
-				   const InVector        &fe_function,
-				   const Function<dim>   &exact_solution,
-				   OutVector             &difference,
-				   const Quadrature<dim> &q,
-				   const NormType        &norm,
-				   const Function<dim>   *weight,
-				   const double           exponent_1)
-{
-				   // we mark the "exponent" parameter
-				   // to this function "const" since
-				   // it is strictly incoming, but we
-				   // need to set it to something
-				   // different later on, if
-				   // necessary, so have a read-write
-				   // version of it:
-  double exponent = exponent_1;
-  
-  const unsigned int        n_q_points   = q.n_quadrature_points;
-  const FiniteElement<dim> &fe           = dof.get_fe();
-  const unsigned int        n_components = fe.n_components();
-  const bool                fe_is_system = (n_components != 1);
-
-  if (weight!=0)
-    {
-      Assert ((weight->n_components==1) || (weight->n_components==n_components),
-	      ExcDimensionMismatch(weight->n_components, n_components));
-    }
-
-  difference.reinit (dof.get_tria().n_active_cells());
-  
-  switch (norm)
-    {
-      case L2_norm:
-      case H1_seminorm:
-      case H1_norm:
-	exponent = 2.;
-	break;
-      case L1_norm:
-	exponent = 1.;
-	break;
-      default:
-	break;
-    }
-  
-  UpdateFlags update_flags = UpdateFlags (update_q_points  |
- 					  update_JxW_values);
-  switch (norm)
-    {
-      case H1_seminorm:
-      case W1p_seminorm:
-      case W1infty_seminorm:
-	update_flags |= UpdateFlags (update_gradients);
-	break;
-      case H1_norm:
-      case W1p_norm:
-      case W1infty_norm:
-	update_flags |= UpdateFlags (update_gradients);
-					 // no break!
-      default:
-	update_flags |= UpdateFlags (update_values);
-	break;
-    }  
-  
-  FEValues<dim> fe_values(mapping, fe, q, update_flags);
-
-  std::vector< Vector<double> >        function_values (n_q_points,
-							Vector<double>(n_components));
-  std::vector<std::vector<Tensor<1,dim> > > function_grads (n_q_points,
-							    std::vector<Tensor<1,dim> >(n_components));
-  std::vector<double> weight_values (n_q_points);
-  std::vector<Vector<double> > weight_vectors (n_q_points, 
-					       Vector<double>(n_components));
-  
-  std::vector<Vector<double> >         psi_values (n_q_points,
-						   Vector<double>(n_components));
-  std::vector<std::vector<Tensor<1,dim> > > psi_grads (n_q_points,
-						       std::vector<Tensor<1,dim> >(n_components));
-  std::vector<double> psi_scalar (n_q_points);
-				   // tmp vector when we use the
-				   // Function<dim> functions for
-				   // scalar functions
-  std::vector<double>         tmp_values (fe_values.n_quadrature_points);
-  std::vector<Tensor<1,dim> > tmp_gradients (fe_values.n_quadrature_points);
-  
- 				   // loop over all cells
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
-						 endc = dof.end();
-  for (unsigned int index=0; cell != endc; ++cell, ++index)
-    {
-      double diff=0;
- 				       // initialize for this cell
-      fe_values.reinit (cell);
-      
-      if (weight!=0)
-	{
-	  if (weight->n_components>1)
-	    weight->vector_value_list (fe_values.get_quadrature_points(),
-				       weight_vectors);
-	  else
-	    {
-	      weight->value_list (fe_values.get_quadrature_points(),
-				  weight_values);
-	      for (unsigned int k=0;k<n_q_points;++k)
-		weight_vectors[k] = weight_values[k];
-	    }
-	} else {
-	  for (unsigned int k=0;k<n_q_points;++k)
-	    weight_vectors[k] = 1.;
-	}
-      
-      
-      if (update_flags & update_values)
-	{
-					   // first compute the exact solution
-					   // (vectors) at the quadrature points
-					   // try to do this as efficient as
-					   // possible by avoiding a second
-					   // virtual function call in case
-					   // the function really has only
-					   // one component
-	  if (fe_is_system)
-	    exact_solution.vector_value_list (fe_values.get_quadrature_points(),
-					      psi_values);
-	  else
-	    {
-	      exact_solution.value_list (fe_values.get_quadrature_points(),
-					 tmp_values);
-	      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);
-	  for (unsigned int q=0; q<n_q_points; ++q)
-	    psi_values[q] -= function_values[q];
-	}
-
-				       // Do the same for gradients, if required
-      if (update_flags & update_gradients)
-	{
-					   // try to be a little clever
-					   // to avoid recursive virtual
-					   // function calls when calling
-					   // @p{gradient_list} for functions
-					   // that are really scalar
-					   // functions
-	  if (fe_is_system)
-	    exact_solution.vector_gradient_list (fe_values.get_quadrature_points(),
-						 psi_grads);
-	  else
-	    {
-	      exact_solution.gradient_list (fe_values.get_quadrature_points(),
-					    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);
-	  for (unsigned int k=0; k<n_components; ++k)
-	    for (unsigned int q=0; q<n_q_points; ++q)
-	      psi_grads[q][k] -= function_grads[q][k];
-	}
-      
-      switch (norm)
-	{
-	  case mean:
-	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
-					     // Compute values in
-					     // quadrature points
-	    for (unsigned int k=0; k<n_components; ++k)
-	      for (unsigned int q=0; q<n_q_points; ++q)
-		  psi_scalar[q] += psi_values[q](k)
-				   * weight_vectors[q](k);
-
-					     // Integrate
-	    diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
-				       fe_values.get_JxW_values().begin(),
-				       0.0);
-	    break;
-	  case Lp_norm:
-	  case L1_norm:
-	  case W1p_norm:
-	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
-					     // Compute values in
-					     // quadrature points
-	    for (unsigned int k=0; k<n_components; ++k)
-	      for (unsigned int q=0; q<n_q_points; ++q)
-		  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(),
-				       0.0);
-					     // Compute the root only,
-					     // if no derivative
-					     // values are added later
-	    if (!(update_flags & update_gradients))
-	      diff = std::pow(diff, 1./exponent);
-	    break;
-	  case L2_norm:
-	  case H1_norm:
-	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
-					     // Compute values in
-					     // quadrature points
-	    for (unsigned int k=0; k<n_components; ++k)
-	      for (unsigned int q=0; q<n_q_points; ++q)
-		  psi_scalar[q] += psi_values[q](k)*psi_values[q](k)
-				   * weight_vectors[q](k);
-
-					     // Integrate
-	    diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
-				       fe_values.get_JxW_values().begin(),
-				       0.0);
-					     // Compute the root only,
-					     // if no derivative
-					     // values are added later
-	    if (norm == L2_norm)
-	      diff=std::sqrt(diff);
-	    break;
-	  case Linfty_norm:
-	  case W1infty_norm:
-	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
-	    for (unsigned int k=0; k<n_components; ++k)
-	      for (unsigned int q=0; q<n_q_points; ++q)
-		{
-		  double newval = std::fabs(psi_values[q](k))
-				  * weight_vectors[q](k);
-		  if (psi_scalar[q]<newval)
-		    psi_scalar[q] = newval;
-		}
-					     // Maximum on one cell
-	    diff = *std::max_element (psi_scalar.begin(), psi_scalar.end());
-	    break;
-	  case H1_seminorm:
-	  case W1p_seminorm:
-	    break;
-	  default:
-	    Assert (false, ExcNotImplemented());
-	    break;
-	}
-
-      switch (norm)
-	{
-	  case W1p_seminorm:
-	  case W1p_norm:
-	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
-	    for (unsigned int k=0; k<n_components; ++k)
-	      for (unsigned int q=0; q<n_q_points; ++q)
-		psi_scalar[q] += std::pow(sqr_point(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);
-	    diff = std::pow(diff, 1./exponent);
-	    break;
-	  case H1_seminorm:
-	  case H1_norm:
-					     // take square of integrand
-	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
-	    for (unsigned int k=0; k<n_components; ++k)
-	      for (unsigned int q=0; q<n_q_points; ++q)
-		psi_scalar[q] += sqr_point(psi_grads[q][k])
-				 * weight_vectors[q](k);
-
- 					     // add seminorm to L_2 norm or
- 					     // to zero
-	    diff += std::inner_product (psi_scalar.begin(), psi_scalar.end(),
-					fe_values.get_JxW_values().begin(),
-					0.0);
-	    diff = std::sqrt(diff);
-	    break;
-	  case W1infty_seminorm:
-	  case W1infty_norm:
-	    Assert(false, ExcNotImplemented());
-	    std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
-	    for (unsigned int k=0; k<n_components; ++k)
-	      for (unsigned int q=0; q<n_q_points; ++q)
-		{
-		  double t = 0.;
-		  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);
-		}
-
-	    for (unsigned int i=0;i<psi_scalar.size();++i)
-	      diff = std::max (diff, psi_scalar[i]);
-	    break;
-	  default:
-	    break;
- 	}
-				       // append result of this cell
- 				       // to the end of the vector
-      difference(index) = diff;
-    }
-}
-
-
-template <int dim, class InVector, class OutVector>
-void
-VectorTools::integrate_difference (const DoFHandler<dim>    &dof,
-				   const InVector           &fe_function,
-				   const Function<dim>      &exact_solution,
-				   OutVector                &difference,
-				   const Quadrature<dim>    &q,
-				   const NormType           &norm,
-				   const Function<dim>      *weight,
-				   const double              exponent)
-{
-  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  static const MappingQ1<dim> mapping;
-  integrate_difference(mapping, dof, fe_function, exact_solution,
-		       difference, q, norm, weight, exponent);
-}
-
-
-
-template <int dim, class InVector>
-void
-VectorTools::point_difference (const DoFHandler<dim> &dof,
-			       const InVector        &fe_function,
-			       const Function<dim>   &exact_function,
-			       Vector<double>        &difference,
-			       const Point<dim>      &point)
-{
-  static const MappingQ1<dim> mapping;
-  const FiniteElement<dim>& fe = dof.get_fe();
-
-  Assert(difference.size() == fe.n_components(),
-	 ExcDimensionMismatch(difference.size(), fe.n_components()));
-
-                                   // first find the cell in which this point
-                                   // is, initialize a quadrature rule with
-                                   // it, and then a FEValues object
-  const typename DoFHandler<dim>::active_cell_iterator
-    cell = GridTools::find_active_cell_around_point (dof, point);
-
-  const Point<dim> unit_point
-    = mapping.transform_real_to_unit_cell(cell, point);
-  Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
-          ExcInternalError());
-
-  const Quadrature<dim> quadrature (std::vector<Point<dim> > (1, unit_point),
-                                    std::vector<double> (1, 1.));
-  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
-  fe_values.reinit(cell);
-
-                                   // then use this to get at the values of
-                                   // the given fe_function at this point
-  std::vector<Vector<double> > u_value(1, Vector<double> (fe.n_components()));
-  fe_values.get_function_values(fe_function, u_value);
-
-  if (fe.n_components() == 1)
-    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);
-}
-
-
-
-template <int dim, class InVector>
-double
-VectorTools::compute_mean_value (const Mapping<dim>    &mapping,
-				 const DoFHandler<dim> &dof,
-				 const Quadrature<dim> &quadrature,
-				 const InVector        &v,
-				 const unsigned int     component)
-{
-  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));
+// explicit instantiations
 
-  typename DoFHandler<dim>::active_cell_iterator c;
-  std::vector<Vector<double> > values(quadrature.n_quadrature_points,
-				      Vector<double> (dof.get_fe().n_components()));
-  
-  double mean = 0.;
-  double area = 0.;
-				   // Compute mean value
-  for (c = dof.begin_active(); c != dof.end(); ++c)
-    {
-      fe.reinit (c);
-      fe.get_function_values(v, values);
-      for (unsigned int k=0; k< quadrature.n_quadrature_points; ++k)
-	{
-	  mean += fe.JxW(k) * values[k](component);
-	  area += fe.JxW(k);
-	};
-    };
-  
-  return (mean/area);
-}
+#define VEC Vector<double>
+#include "vectors.instance.h"
+#undef VEC
 
+#define VEC Vector<float>
+#include "vectors.instance.h"
+#undef VEC
 
-template <int dim, class InVector>
-double
-VectorTools::compute_mean_value (const DoFHandler<dim> &dof,
-				 const Quadrature<dim> &quadrature,
-				 const InVector        &v,
-				 const unsigned int     component)
-{
-  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
-  static const MappingQ1<dim> mapping;
-  return compute_mean_value(mapping, dof, quadrature, v, component);
-}
+#define VEC BlockVector<double>
+#include "vectors.instance.h"
+#undef VEC
 
-
-// explicit instantiations
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const Mapping<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
- Vector<double>&);
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
- Vector<double>&);
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const Mapping<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
- Vector<float>&);
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
- Vector<float>&);
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const Mapping<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
- BlockVector<double>&);
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
- BlockVector<double>&);
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const Mapping<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
- BlockVector<float>&);
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
- BlockVector<float>&);
-
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const FullMatrix<double>&,
- const Vector<double>&,
- Vector<double>&);
-template
-void VectorTools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const FullMatrix<double>&,
- const BlockVector<double>&,
- BlockVector<double>&);
+#define VEC BlockVector<float>
+#include "vectors.instance.h"
+#undef VEC
 
 template
 void VectorTools::project<deal_II_dimension>
@@ -1721,6 +41,7 @@ void VectorTools::project<deal_II_dimension>
  const bool,
  const Quadrature<deal_II_dimension-1> &,
  const bool);
+
 template
 void VectorTools::create_right_hand_side<deal_II_dimension>
 (const Mapping<deal_II_dimension>    &,
@@ -1762,204 +83,6 @@ void VectorTools::interpolate_boundary_values<deal_II_dimension> (
   const Function<deal_II_dimension>   &,
   std::map<unsigned int,double>       &,
   const std::vector<bool>    &);
-template
-void VectorTools::integrate_difference<deal_II_dimension> (
-  const Mapping<deal_II_dimension>    &,
-  const DoFHandler<deal_II_dimension> &,
-  const Vector<double>                &,
-  const Function<deal_II_dimension>   &,
-  Vector<float>                       &,
-  const Quadrature<deal_II_dimension> &,
-  const NormType                      &,
-  const Function<deal_II_dimension>   *,
-  const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &,
- const Vector<double>                &,
- const Function<deal_II_dimension>   &,
- Vector<float>                       &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const Mapping<deal_II_dimension>    &,
- const DoFHandler<deal_II_dimension> &,
- const Vector<double>                &,
- const Function<deal_II_dimension>   &,
- Vector<double>                      &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &,
- const Vector<double>                &,
- const Function<deal_II_dimension>   &,
- Vector<double>                      &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const Mapping<deal_II_dimension>    &,
- const DoFHandler<deal_II_dimension> &,
- const Vector<float>                 &,
- const Function<deal_II_dimension>   &,
- Vector<float>                       &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &,
- const Vector<float>                 &,
- const Function<deal_II_dimension>   &,
- Vector<float>                       &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const Mapping<deal_II_dimension>    &,
- const DoFHandler<deal_II_dimension> &,
- const Vector<float>                 &,
- const Function<deal_II_dimension>   &,
- Vector<double>                      &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &,
- const Vector<float>                 &,
- const Function<deal_II_dimension>   &,
- Vector<double>                      &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const Mapping<deal_II_dimension>    &,
- const DoFHandler<deal_II_dimension> &,
- const BlockVector<double>                &,
- const Function<deal_II_dimension>   &,
- Vector<float>                       &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &,
- const BlockVector<double>                &,
- const Function<deal_II_dimension>   &,
- Vector<float>                       &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const Mapping<deal_II_dimension>    &,
- const DoFHandler<deal_II_dimension> &,
- const BlockVector<double>                &,
- const Function<deal_II_dimension>   &,
- Vector<double>                      &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &,
- const BlockVector<double>                &,
- const Function<deal_II_dimension>   &,
- Vector<double>                      &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const Mapping<deal_II_dimension>    &,
- const DoFHandler<deal_II_dimension> &,
- const BlockVector<float>                 &,
- const Function<deal_II_dimension>   &,
- Vector<float>                       &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &,
- const BlockVector<float>                 &,
- const Function<deal_II_dimension>   &,
- Vector<float>                       &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const Mapping<deal_II_dimension>    &,
- const DoFHandler<deal_II_dimension> &,
- const BlockVector<float>                 &,
- const Function<deal_II_dimension>   &,
- Vector<double>                      &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-template
-void VectorTools::integrate_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &,
- const BlockVector<float>                 &,
- const Function<deal_II_dimension>   &,
- Vector<double>                      &,
- const Quadrature<deal_II_dimension> &,
- const NormType                      &,
- const Function<deal_II_dimension>   *,
- const double);
-
-template
-void VectorTools::point_difference<deal_II_dimension> (
-  const DoFHandler<deal_II_dimension>&,
-  const Vector<double>&,
-  const Function<deal_II_dimension>&,
-  Vector<double>&,
-  const Point<deal_II_dimension>&);
-template
-void VectorTools::point_difference<deal_II_dimension> (
-  const DoFHandler<deal_II_dimension>&,
-  const Vector<float>&,
-  const Function<deal_II_dimension>&,
-  Vector<double>&,
-  const Point<deal_II_dimension>&);
-template
-void VectorTools::point_difference<deal_II_dimension> (
-  const DoFHandler<deal_II_dimension>&,
-  const BlockVector<double>&,
-  const Function<deal_II_dimension>&,
-  Vector<double>&,
-  const Point<deal_II_dimension>&);
-template
-void VectorTools::point_difference<deal_II_dimension> (
-  const DoFHandler<deal_II_dimension>&,
-  const BlockVector<float>&,
-  const Function<deal_II_dimension>&,
-  Vector<double>&,
-  const Point<deal_II_dimension>&);
-
 
 #if deal_II_dimension != 1
 template
@@ -1978,63 +101,6 @@ void VectorTools::project_boundary_values<deal_II_dimension>
  const Quadrature<deal_II_dimension-1>&,
  std::map<unsigned int,double>        &);
 
-template
-double VectorTools::compute_mean_value<deal_II_dimension>
-(const Mapping<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const Quadrature<deal_II_dimension>&,
- const Vector<double>&,
- const unsigned int);
-template
-double VectorTools::compute_mean_value<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const Quadrature<deal_II_dimension>&,
- const Vector<double>&,
- const unsigned int);
-
-template
-double VectorTools::compute_mean_value<deal_II_dimension>
-(const Mapping<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const Quadrature<deal_II_dimension>&,
- const BlockVector<double>&,
- const unsigned int);
-template
-double VectorTools::compute_mean_value<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const Quadrature<deal_II_dimension>&,
- const BlockVector<double>&,
- const unsigned int);
-
-template
-double VectorTools::compute_mean_value<deal_II_dimension>
-(const Mapping<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const Quadrature<deal_II_dimension>&,
- const Vector<float>&,
- const unsigned int);
-template
-double VectorTools::compute_mean_value<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const Quadrature<deal_II_dimension>&,
- const Vector<float>&,
- const unsigned int);
-
-template
-double VectorTools::compute_mean_value<deal_II_dimension>
-(const Mapping<deal_II_dimension>&,
- const DoFHandler<deal_II_dimension>&,
- const Quadrature<deal_II_dimension>&,
- const BlockVector<float>&,
- const unsigned int);
-template
-double VectorTools::compute_mean_value<deal_II_dimension>
-(const DoFHandler<deal_II_dimension>&,
- const Quadrature<deal_II_dimension>&,
- const BlockVector<float>&,
- const unsigned int);
-
-
 
 // the following two functions are not derived from a template in 1d
 // and thus need no explicit instantiation
diff --git a/deal.II/deal.II/source/numerics/vectors.instance.h b/deal.II/deal.II/source/numerics/vectors.instance.h
new file mode 100644
index 0000000000..a2d1464b22
--- /dev/null
+++ b/deal.II/deal.II/source/numerics/vectors.instance.h
@@ -0,0 +1,103 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+// Instantiations of functions in vectors.cc
+
+// vectors.cc defines the vector type VEC
+
+template
+void VectorTools::interpolate<deal_II_dimension>
+(const Mapping<deal_II_dimension>&,
+ const DoFHandler<deal_II_dimension>&,
+ const Function<deal_II_dimension>&,
+ VEC&);
+template
+void VectorTools::interpolate<deal_II_dimension>
+(const DoFHandler<deal_II_dimension>&,
+ const Function<deal_II_dimension>&,
+ VEC&);
+
+// Should these be instantiated for every combination of two types?
+template
+void VectorTools::interpolate<deal_II_dimension>
+(const DoFHandler<deal_II_dimension>&,
+ const DoFHandler<deal_II_dimension>&,
+ const FullMatrix<double>&,
+ const VEC&,
+ VEC&);
+
+
+template
+void VectorTools::integrate_difference<deal_II_dimension>
+(const DoFHandler<deal_II_dimension>&,
+ const VEC&,
+ const Function<deal_II_dimension>&,
+ Vector<double>&,
+ const Quadrature<deal_II_dimension>&,
+ const NormType&,
+ const Function<deal_II_dimension>*,
+ const double);
+template
+void VectorTools::integrate_difference<deal_II_dimension>
+(const DoFHandler<deal_II_dimension>&,
+ const VEC&,
+ const Function<deal_II_dimension>&,
+ Vector<float>&,
+ const Quadrature<deal_II_dimension>&,
+ const NormType&,
+ const Function<deal_II_dimension>*,
+ const double);
+template
+void VectorTools::integrate_difference<deal_II_dimension>
+(const Mapping<deal_II_dimension>&,
+ const DoFHandler<deal_II_dimension>&,
+ const VEC&,
+ const Function<deal_II_dimension>&,
+ Vector<double>&,
+ const Quadrature<deal_II_dimension>&,
+ const NormType&,
+ const Function<deal_II_dimension>*,
+ const double);
+template
+void VectorTools::integrate_difference<deal_II_dimension>
+(const Mapping<deal_II_dimension>&,
+ const DoFHandler<deal_II_dimension>&,
+ const VEC&,
+ const Function<deal_II_dimension>&,
+ Vector<float>&,
+ const Quadrature<deal_II_dimension>&,
+ const NormType&,
+ const Function<deal_II_dimension>*,
+ const double);
+
+template
+void VectorTools::point_difference<deal_II_dimension> (
+  const DoFHandler<deal_II_dimension>&,
+  const VEC&,
+  const Function<deal_II_dimension>&,
+  Vector<double>&,
+  const Point<deal_II_dimension>&);
+
+template
+double VectorTools::compute_mean_value<deal_II_dimension>
+(const Mapping<deal_II_dimension>&,
+ const DoFHandler<deal_II_dimension>&,
+ const Quadrature<deal_II_dimension>&,
+ const VEC&,
+ const unsigned int);
+template
+double VectorTools::compute_mean_value<deal_II_dimension>
+(const DoFHandler<deal_II_dimension>&,
+ const Quadrature<deal_II_dimension>&,
+ const VEC&,
+ const unsigned int);