]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
move instatiations to vectors.templates.h; thanks to Dima Sorkin
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 30 Aug 2005 09:48:05 +0000 (09:48 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 30 Aug 2005 09:48:05 +0000 (09:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@11339 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/include/numerics/vectors.templates.h [new file with mode: 0644]
deal.II/deal.II/source/numerics/vectors.cc
deal.II/deal.II/source/numerics/vectors.instance.h [new file with mode: 0644]

index 22447c2012d5bb2ae0c5bbf22c4f3f9d51535412..2b37acf61a209d68be23a49cfaea08bbdf7bb707 100644 (file)
@@ -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 (file)
index 0000000..4f07b68
--- /dev/null
@@ -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
index 2bbad617ac37e6b54752c3aff7abb1a994504dfa..3dbc4d177932c1a6e72afd33b9948825808ae734 100644 (file)
 //
 //---------------------------------------------------------------------------
 
+#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 (file)
index 0000000..a2d1464
--- /dev/null
@@ -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);

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.