From 6e8b912c4ed34e33b0d7346bc6649423b8399e57 Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 29 Aug 2007 21:46:46 +0000
Subject: [PATCH] Update comment and indentation style

git-svn-id: https://svn.dealii.org/trunk@15092 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/base/include/base/parsed_function.h   | 387 +++++++------
 deal.II/base/source/parsed_function.cc        | 113 ++--
 .../include/numerics/fe_field_function.h      | 508 +++++++++++-------
 .../numerics/fe_field_function.templates.h    | 338 +++++++-----
 .../source/numerics/fe_field_function.cc      |  23 +-
 5 files changed, 807 insertions(+), 562 deletions(-)

diff --git a/deal.II/base/include/base/parsed_function.h b/deal.II/base/include/base/parsed_function.h
index 69da46e8dc..913a5db348 100644
--- a/deal.II/base/include/base/parsed_function.h
+++ b/deal.II/base/include/base/parsed_function.h
@@ -20,178 +20,241 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace Functions {
-  /** Friendly interface to the FunctionParser class. This class is
-      meant as a wrapper for the FunctionParser class. It provides two
-      methods to declare and parse a ParameterHandler object and creates
-      the Function object declared in the parameter file. This class is
-      derived from the AutoDerivativeFunction class, so you don't need
-      to specify derivatives. An example of usage of this class is as follows:
-    
-      \code
-      // A parameter handler
-      ParameterHandler prm;
-    
-      // Declare a section for the function we need
-      prm.enter_subsection("My vector function");
-      ParsedFunction<dim>::declare_parameters(prm, dim);
-      prm.leave_subsection();
-    
-      // Create a ParsedFunction
-      ParsedFunction<dim> my_vector_function(dim);
-    
-      // Parse an input file.
-      prm.read_input(some_input_file);
-    
-      // Initialize the ParsedFunction object with the given file
-      prm.enter_subsection("My vector function");
-      my_vector_function.parse_parameters(prm);
-      prm.leave_subsection();
-
-      \endcode
-    
-      And here is an example of how the input parameter could look like
-      (see the documentation of the FunctionParser class for a detailed
-      description of the syntax of the function definition):
-    
-      \code
-    
-      # A test two dimensional vector function, depending on time
-      subsection My vector function
-      set Function constants  = kappa=.1, lambda=2.
-      set Function expression = if(y>.5, kappa*x*(1-x),0); t^2*cos(lambda*pi*x)
-      set Variable names      = x,y,t
-      end
-    
-      \endcode
-    
-      \ingroup functions
-      \author Luca Heltai, 2006
-  */
+namespace Functions
+{
+/**
+ *   Friendly interface to the FunctionParser class. This class is
+ *   meant as a wrapper for the FunctionParser class. It provides two
+ *   methods to declare and parse a ParameterHandler object and creates
+ *   the Function object declared in the parameter file. This class is
+ *   derived from the AutoDerivativeFunction class, so you don't need
+ *   to specify derivatives. An example of usage of this class is as follows:
+ * 
+ *   @code
+ *   // A parameter handler
+ *   ParameterHandler prm;
+ * 
+ *   // Declare a section for the function we need
+ *   prm.enter_subsection("My vector function");
+ *   ParsedFunction<dim>::declare_parameters(prm, dim);
+ *   prm.leave_subsection();
+ * 
+ *   // Create a ParsedFunction
+ *   ParsedFunction<dim> my_vector_function(dim);
+ * 
+ *   // Parse an input file.
+ *   prm.read_input(some_input_file);
+ * 
+ *   // Initialize the ParsedFunction object with the given file
+ *   prm.enter_subsection("My vector function");
+ *   my_vector_function.parse_parameters(prm);
+ *   prm.leave_subsection();
+ *
+ *   @endcode
+ * 
+ *   And here is an example of how the input parameter could look like
+ *   (see the documentation of the FunctionParser class for a detailed
+ *   description of the syntax of the function definition):
+ * 
+ *   @code
+ * 
+ *   # A test two dimensional vector function, depending on time
+ *   subsection My vector function
+ *   set Function constants  = kappa=.1, lambda=2.
+ *   set Function expression = if(y>.5, kappa*x*(1-x),0); t^2*cos(lambda*pi*x)
+ *   set Variable names      = x,y,t
+ *   end
+ * 
+ *   @endcode
+ * 
+ *   @ingroup functions
+ *   @author Luca Heltai, 2006
+ */
   template <int dim>
   class ParsedFunction :  public AutoDerivativeFunction<dim> 
   {
-  public:
-    /** Construct a vector function. The vector function which is
-	generated has @p n_components components (defaults to 1). The parameter
-	@p h is used to initialize the AutoDerivativeFunction class from
-	which this class is derived. */
-    ParsedFunction (const unsigned int n_components = 1, const double h=1e-8);
+    public:
+				       /**
+					* Construct a vector
+					* function. The vector
+					* function which is generated
+					* has @p n_components
+					* components (defaults to
+					* 1). The parameter @p h is
+					* used to initialize the
+					* AutoDerivativeFunction class
+					* from which this class is
+					* derived.
+					*/
+      ParsedFunction (const unsigned int n_components = 1, const double h=1e-8);
   
-    /** Declare parameters needed by this class. The additional
-	parameter @p n_components is used to generate the right code according
-	to the number of components of the function that will parse this
-	ParameterHandler. If the number of components which is parsed
-	does not match the number of components of this object, an
-	assertion is thrown and the program is aborted. 
-      
-	The default behavior for this class is to declare the following
-	entries:
-      
-	\code
-      
-	set Function constants  = 
-	set Function expression = 0
-	set Variable names      = x,y,t
-      
-	\endcode
-      
-    */
-    static void declare_parameters(ParameterHandler &prm, 
-				   const unsigned int n_components = 1);
+				       /**
+					* Declare parameters needed by
+					* this class. The additional
+					* parameter @p n_components is
+					* used to generate the right
+					* code according to the number
+					* of components of the
+					* function that will parse
+					* this ParameterHandler. If
+					* the number of components
+					* which is parsed does not
+					* match the number of
+					* components of this object,
+					* an assertion is thrown and
+					* the program is aborted.  The
+					* default behavior for this
+					* class is to declare the
+					* following entries:
+					*
+					*  @code
+					*
+					*  set Function constants  = 
+					*  set Function expression = 0
+					*  set Variable names      = x,y,t
+					*
+					*  @endcode
+					*
+					*/
+      static void declare_parameters(ParameterHandler &prm, 
+				     const unsigned int n_components = 1);
   
-    /** Parse parameters needed by this class.  If the number of
-	components which is parsed does not match the number of
-	components of this object, an assertion is thrown and the
-	program is aborted. 
-      
-	In order for the class to function properly, we follow the same
-	convenctions declared in the FunctionParser class (look there
-	for a detailed description of the syntax for function
-	declarations).
-      
-	The three variables that can be parsed from a parameter file are
-	the following:
-      
-	\code
-      
-	set Function constants  = 
-	set Function expression = 
-	set Variable names      = 
-      
-	\endcode
-
-	Function constants is a collection of pairs in the form
-	name=value, separated by commas, for example:
-
-	\code
-      
-	set Function constants = lambda=1. , alpha=2., gamma=3.
-      
-	\endcode
-      
-	These constants can be used in the declaration of the function
-	expression, which follows the convention of the FunctionParser
-	class. In order to specify vector functions, semicolons have to
-	be used to separate the different components, e.g.:
-      
-	\code
-      
-	set Function expression = cos(pi*x) ; cos(pi*y) 
-
-	\endcode
-      
-	The variable names entry can be used to customize the name of
-	the variables used in the Function. It defaults to
-
-	\code
-      
-	set Variable names      = x,t
-      
-	\endcode
-      
-	for one dimensional problems, 
-
-	\code
-      
-	set Variable names      = x,y,t
-      
-	\endcode
-      
-	for two dimensional problems and
-
-	\code
-      
-	set Variable names      = x,y,z,t
-      
-	\endcode
-      
-	for three dimensional problems.
-      
-	The time variable can be set according to specifications in the
-	FunctionTime class.
-      
-    */
-    void parse_parameters(ParameterHandler &prm);
+				       /**
+					* Parse parameters needed by
+					* this class.  If the number
+					* of components which is
+					* parsed does not match the
+					* number of components of this
+					* object, an assertion is
+					* thrown and the program is
+					* aborted.  In order for the
+					* class to function properly,
+					* we follow the same
+					* convenctions declared in the
+					* FunctionParser class (look
+					* there for a detailed
+					* description of the syntax
+					* for function declarations).
+					*
+					* The three variables that can
+					* be parsed from a parameter
+					* file are the following:
+					*
+					*  @code
+					*
+					*  set Function constants  = 
+					*  set Function expression = 
+					*  set Variable names      = 
+					*
+					*  @endcode
+					*
+					*  Function constants is a
+					*  collection of pairs in the
+					*  form name=value, separated
+					*  by commas, for example:
+					*
+					*  @code
+					*
+					*  set Function constants = lambda=1. , alpha=2., gamma=3.
+					*
+					*  @endcode
+					*
+					*  These constants can be used
+					*  in the declaration of the
+					*  function expression, which
+					*  follows the convention of
+					*  the FunctionParser
+					*  class. In order to specify
+					*  vector functions,
+					*  semicolons have to be used
+					*  to separate the different
+					*  components, e.g.:
+					*
+					*  @code
+					*
+					*  set Function expression = cos(pi*x) ; cos(pi*y) 
+					*
+					*  @endcode
+					*
+					*  The variable names entry
+					*  can be used to customize
+					*  the name of the variables
+					*  used in the Function. It
+					*  defaults to
+					*
+					*  @code
+					*
+					*  set Variable names      = x,t
+					*
+					*  @endcode
+					*
+					*  for one dimensional problems, 
+					*
+					*  @code
+					*
+					*  set Variable names      = x,y,t
+					*
+					*  @endcode
+					*
+					*  for two dimensional problems and
+					*
+					*  @code
+					*
+					*  set Variable names      = x,y,z,t
+					*
+					*  @endcode
+					*
+					*  for three dimensional problems.
+					*
+					*  The time variable can be
+					*  set according to
+					*  specifications in the
+					*  FunctionTime base class.
+					*/
+      void parse_parameters(ParameterHandler &prm);
   
-    /** Get one value at the given point. */ 
-    virtual void vector_value (const Point<dim> &p,
-			       Vector<double>   &values) const;
+				       /**
+					* Get one value at the given
+					* point.
+					*/ 
+      virtual void vector_value (const Point<dim> &p,
+				 Vector<double>   &values) const;
   
-    /** Return the value of the function at the given point. Unless
-	there is only one component (i.e. the function is scalar), you
-	should state the component you want to have evaluated; it
-	defaults to zero, i.e. the first component. */
-    virtual double value (const Point< dim > &    p,
-			  const unsigned int  component = 0)    const;
+				       /**
+					* Return the value of the
+					* function at the given
+					* point. Unless there is only
+					* one component (i.e. the
+					* function is scalar), you
+					* should state the component
+					* you want to have evaluated;
+					* it defaults to zero,
+					* i.e. the first component.
+					*/
+      virtual double value (const Point< dim > &    p,
+			    const unsigned int  component = 0)    const;
   
-    /** We need to overwrite this to set the time also in the accessor
-	FunctionParser<dim>. */
-    virtual void set_time(const double newtime);
-  private:
-    FunctionParser<dim> function_object;
+				       /**
+					* Set the time to a specific
+					* value for time-dependent
+					* functions.
+					*
+					* We need to overwrite this to
+					* set the time also in the
+					* accessor
+					* FunctionParser<dim>.
+					*/
+      virtual void set_time(const double newtime);
+      
+    private:
+				       /**
+					* The object with which we do
+					* computations.
+					*/
+      FunctionParser<dim> function_object;
   };
 }
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
diff --git a/deal.II/base/source/parsed_function.cc b/deal.II/base/source/parsed_function.cc
index 995eb64bb6..3ff828199c 100644
--- a/deal.II/base/source/parsed_function.cc
+++ b/deal.II/base/source/parsed_function.cc
@@ -17,45 +17,56 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace Functions {
   template <int dim>
-  ParsedFunction<dim>::ParsedFunction (const unsigned int n_components, const double h) : 
-    AutoDerivativeFunction<dim>(h, n_components),
-    function_object(n_components)
-  {
-  }
+  ParsedFunction<dim>::ParsedFunction (const unsigned int n_components, const double h)
+		  : 
+		  AutoDerivativeFunction<dim>(h, n_components),
+		  function_object(n_components)
+  {}
 
+
+  
   template <int dim>
-  void ParsedFunction<dim>::declare_parameters(ParameterHandler &prm, unsigned int n_components) 
+  void
+  ParsedFunction<dim>::declare_parameters(ParameterHandler  &prm,
+					  const unsigned int n_components) 
   {
     Assert(n_components > 0, ExcZero());
 
     std::string vnames;
-    switch (dim) {
-    case 1:
-      vnames = "x,t";
-      break;
-    case 2:
-      vnames = "x,y,t";
-      break;
-    case 3:
-      vnames = "x,y,z,t";
-      break;
-    default:
-      AssertThrow(false, ExcNotImplemented());
-      break;
-    }
+    switch (dim)
+      {
+	case 1:
+	      vnames = "x,t";
+	      break;
+	case 2:
+	      vnames = "x,y,t";
+	      break;
+	case 3:
+	      vnames = "x,y,z,t";
+	      break;
+	default:
+	      AssertThrow(false, ExcNotImplemented());
+	      break;
+      }
     prm.declare_entry("Variable names", vnames, Patterns::Anything(), 
-		      "The name of the variables as they will be used in the function, separated by ','.");
-    // The expression of the function
+		      "The name of the variables as they will be used in the "
+		      "function, separated by ','.");
+
+				     // The expression of the function
     std::string expr = "0";
-    for(unsigned int i=1; i<n_components; ++i)
+    for (unsigned int i=1; i<n_components; ++i)
       expr += "; 0";
 
     prm.declare_entry("Function expression", expr, Patterns::Anything(),
-		      "Separate vector valued expressions by ';' as ',' is used internally by the function parser.");
+		      "Separate vector valued expressions by ';' as ',' "
+		      "is used internally by the function parser.");
     prm.declare_entry("Function constants", "", Patterns::Anything(),
-		      "Any constant used inside the function which is not a variable name.");
+		      "Any constant used inside the function which is not "
+		      "a variable name.");
   }
 
+
+  
   template <int dim>
   void ParsedFunction<dim>::parse_parameters(ParameterHandler &prm) 
   {
@@ -66,33 +77,38 @@ namespace Functions {
     std::vector<std::string> const_list = 
       Utilities::split_string_list(constants_list, ',');
     std::map<std::string, double> constants;
-    for(unsigned int i = 0; i < const_list.size(); ++i) {
-      std::vector<std::string> this_c = 
-	Utilities::split_string_list(const_list[i], '=');
-      AssertThrow(this_c.size() == 2, ExcMessage("Invalid format"));
-      double tmp;
-      AssertThrow( sscanf(this_c[1].c_str(), "%lf", &tmp), ExcMessage("Double number?"));
-      constants[this_c[0]] = tmp;
-    }
+    for(unsigned int i = 0; i < const_list.size(); ++i)
+      {
+	std::vector<std::string> this_c = 
+	  Utilities::split_string_list(const_list[i], '=');
+	AssertThrow(this_c.size() == 2, ExcMessage("Invalid format"));
+	double tmp;
+	AssertThrow( sscanf(this_c[1].c_str(), "%lf", &tmp),
+		     ExcMessage("Double number?"));
+	constants[this_c[0]] = tmp;
+      }
     
     constants["pi"] = deal_II_numbers::PI;
     constants["Pi"] = deal_II_numbers::PI;
 
     unsigned int nn = (Utilities::split_string_list(vnames)).size();
-    switch (nn) {
-    case dim:
-      // Time independent function
-      function_object.initialize(vnames, expression, constants); 
-      break;
-    case dim+1:
-      // Time dependent function
-      function_object.initialize(vnames, expression, constants, true);
-      break;
-    default:
-      AssertThrow(false, ExcMessage("Not the correct size. Check your code."));
-    }
+    switch (nn)
+      {
+	case dim:
+					       // Time independent function
+	      function_object.initialize(vnames, expression, constants); 
+	      break;
+	case dim+1:
+					       // Time dependent function
+	      function_object.initialize(vnames, expression, constants, true);
+	      break;
+	default:
+	      AssertThrow(false, ExcMessage("Not the correct size. Check your code."));
+      }
   }
 
+
+  
   template <int dim>
   void ParsedFunction<dim>::vector_value (const Point<dim> &p,
 					  Vector<double>   &values) const 
@@ -100,6 +116,8 @@ namespace Functions {
     function_object.vector_value(p, values);
   }
 
+
+  
   template <int dim>
   double ParsedFunction<dim>::value (const Point<dim>   &p,
 				     unsigned int comp) const
@@ -107,6 +125,8 @@ namespace Functions {
     return function_object.value(p, comp);
   }
 
+
+  
   template <int dim>
   void ParsedFunction<dim>::set_time (const double newtime)
   { 
@@ -114,7 +134,8 @@ namespace Functions {
     AutoDerivativeFunction<dim>::set_time(newtime);
   }
 
-  // Explicit instantiations
+  
+// Explicit instantiations
   template class ParsedFunction<1>;
   template class ParsedFunction<2>;
   template class ParsedFunction<3>;
diff --git a/deal.II/deal.II/include/numerics/fe_field_function.h b/deal.II/deal.II/include/numerics/fe_field_function.h
index 3137e18e31..61818e4fd3 100644
--- a/deal.II/deal.II/include/numerics/fe_field_function.h
+++ b/deal.II/deal.II/include/numerics/fe_field_function.h
@@ -25,226 +25,342 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace Functions {
+namespace Functions
+{
 
-  /** This is an interpolation function for the given dof handler and
-      the given solution vector. The points at which this function can
-      be evaluated MUST be inside the domain of the dof handler, but
-      except from this, no other requirement is given. This function is
-      rather slow, as it needs to construct a quadrature object for the
-      point (or set of points) where you want to evaluate your finite
-      element function. In order to do so, it needs to find out where
-      the points lie.
-    
-      If you know in advance in which cell your points lye, you can
-      accelerate things a bit, by calling set_active_cell before
-      asking for values or gradients of the function. If you don't do
-      this, and your points don't lie in the cell that is currently
-      stored, the function GridTools::find_cell_around_point is called
-      to find out where the point is. You can specify an optional
-      mapping to use when looking for points in the grid. If you don't
-      do so, this function uses a Q1 mapping.
-    
-      Once the FEFieldFunction knows where the points lie, it creates a
-      quadrature formula for those points, and calls
-      FEValues::get_function_values or FEValues::get_function_grads with
-      the given quadrature points.
-
-      If you only need the quadrature points but not the values of the
-      finite element function (you might want this for the adjoint
-      interpolation), you can also use the function @p
-      compute_point_locations alone.
-    
-      An example of how to use this function is the following:
-    
-      \code
-    
-      // Generate two triangulations
-      Triangulation<dim> tria_1;
-      Triangulation<dim> tria_2;
-    
-      // Read the triangulations from files, or build them up, or get
-      // them from some place...  Assume that tria_2 is *entirely*
-      // included in tria_1
-      ...
-
-      // Associate a dofhandler and a solution to the first
-      // triangulation
-      DoFHandler<dim> dh1(tria_1);
-      Vector<double> solution_1;
-    
-      // Do the same with the second
-      DoFHandler<dim> dh2;
-      Vector<double> solution_2;
-    
-      // Setup the system, assemble matrices, solve problems and get the
-      // nobel prize on the first domain...
-      ...
-    
-      // Now project it to the second domain
-      FEFieldFunction<dim> fe_function_1 (dh_1, solution_1);
-      VectorTools::project(dh_2, constraints_2, quad, fe_function_1, solution_2);
-    
-      // Or interpolate it...
-      Vector<double> solution_3;
-      VectorTools::interpolate(dh_2, fe_function_1, solution_3);
-    
-      \endcode
-    
-      The snippet of code above will work assuming that the second
-      triangulation is entirely included in the first one. 
-    
-      FEFieldFunction is designed to be an easy way to get the results of
-      your computations across different, possibly non matching,
-      grids. No knowledge of the location of the points is assumed in
-      this class, which makes it rely entirely on the
-      GridTools::find_active_cell_around_point utility for its
-      job. However the class can be fed an "educated guess" of where the
-      points that will be computed actually are by using the
-      FEFieldFunction::set_active_cell method, so if you have a smart way to
-      tell where your points are, you will save a lot of computational
-      time by letting this class know.
-
-      An optimization based on a caching mechanism was used by the
-      author of this class for the implementation of a Finite Element
-      Immersed Boundary Method.
-    
-      \addtogroup functions
-
-      \author Luca Heltai, 2006    
-      
-      \todo Add hp functionality
-  */
+/**
+ * This is an interpolation function for the given dof handler and
+ * the given solution vector. The points at which this function can
+ * be evaluated MUST be inside the domain of the dof handler, but
+ * except from this, no other requirement is given. This function is
+ * rather slow, as it needs to construct a quadrature object for the
+ * point (or set of points) where you want to evaluate your finite
+ * element function. In order to do so, it needs to find out where
+ * the points lie.
+ *
+ * If you know in advance in which cell your points lye, you can
+ * accelerate things a bit, by calling set_active_cell before
+ * asking for values or gradients of the function. If you don't do
+ * this, and your points don't lie in the cell that is currently
+ * stored, the function GridTools::find_cell_around_point is called
+ * to find out where the point is. You can specify an optional
+ * mapping to use when looking for points in the grid. If you don't
+ * do so, this function uses a Q1 mapping.
+ *
+ * Once the FEFieldFunction knows where the points lie, it creates a
+ * quadrature formula for those points, and calls
+ * FEValues::get_function_values or FEValues::get_function_grads with
+ * the given quadrature points.
+ *
+ * If you only need the quadrature points but not the values of the
+ * finite element function (you might want this for the adjoint
+ * interpolation), you can also use the function @p
+ * compute_point_locations alone.
+ *
+ * An example of how to use this function is the following:
+ *    
+ * \code
+ * 
+ * // Generate two triangulations
+ * Triangulation<dim> tria_1;
+ * Triangulation<dim> tria_2;
+ *
+ * // Read the triangulations from files, or build them up, or get
+ * // them from some place...  Assume that tria_2 is *entirely*
+ * // included in tria_1
+ * ...
+ *
+ * // Associate a dofhandler and a solution to the first
+ * // triangulation
+ * DoFHandler<dim> dh1(tria_1);
+ * Vector<double> solution_1;
+ *
+ * // Do the same with the second
+ * DoFHandler<dim> dh2;
+ * Vector<double> solution_2;
+ *
+ * // Setup the system, assemble matrices, solve problems and get the
+ * // nobel prize on the first domain...
+ * ...
+ *
+ * // Now project it to the second domain
+ * FEFieldFunction<dim> fe_function_1 (dh_1, solution_1);
+ * VectorTools::project(dh_2, constraints_2, quad, fe_function_1, solution_2);
+ *
+ * // Or interpolate it...
+ * Vector<double> solution_3;
+ * VectorTools::interpolate(dh_2, fe_function_1, solution_3);
+ *
+ * \endcode
+ *
+ * The snippet of code above will work assuming that the second
+ * triangulation is entirely included in the first one. 
+ *
+ * FEFieldFunction is designed to be an easy way to get the results of
+ * your computations across different, possibly non matching,
+ * grids. No knowledge of the location of the points is assumed in
+ * this class, which makes it rely entirely on the
+ * GridTools::find_active_cell_around_point utility for its
+ * job. However the class can be fed an "educated guess" of where the
+ * points that will be computed actually are by using the
+ * FEFieldFunction::set_active_cell method, so if you have a smart way to
+ * tell where your points are, you will save a lot of computational
+ * time by letting this class know.
+ *
+ * An optimization based on a caching mechanism was used by the
+ * author of this class for the implementation of a Finite Element
+ * Immersed Boundary Method.
+ *
+ *  \addtogroup functions
+ *
+ *  \author Luca Heltai, 2006    
+ *  
+ *  \todo Add hp functionality
+ */
   template <int dim, 
 	    typename DH=DoFHandler<dim>,
 	    typename VECTOR=Vector<double> >
   class FEFieldFunction :  public Function<dim> 
   {
-  public:
-    /** Construct a vector function. A smart pointers is stored to the
-	dof handler, so you have to make sure that it make sense for
-	the entire lifetime of this object. The number of components
-	of this functions is equal to the number of components of the
-	finite element object. If a mapping is specified, that is what
-	is used to find out where the points lay. Otherwise the
-	standard Q1 mapping is used. */
-    FEFieldFunction (const DH &dh, const VECTOR &data_vector, 
-		     const Mapping<dim> &mapping = StaticMappingQ1<dim>::mapping);
+    public:
+				       /**
+					* Construct a vector
+					* function. A smart pointers
+					* is stored to the dof
+					* handler, so you have to make
+					* sure that it make sense for
+					* the entire lifetime of this
+					* object. The number of
+					* components of this functions
+					* is equal to the number of
+					* components of the finite
+					* element object. If a mapping
+					* is specified, that is what
+					* is used to find out where
+					* the points lay. Otherwise
+					* the standard Q1 mapping is
+					* used.
+					*/
+      FEFieldFunction (const DH           &dh,
+		       const VECTOR       &data_vector, 
+		       const Mapping<dim> &mapping = StaticMappingQ1<dim>::mapping);
   
-    /** Set the current cell. If you know in advance where your points
-	lie, you can tell this object by calling this function. This
-	will speed things up a little. */
-    inline void set_active_cell(typename DH::active_cell_iterator &newcell);
+				       /**
+					* Set the current cell. If you
+					* know in advance where your
+					* points lie, you can tell
+					* this object by calling this
+					* function. This will speed
+					* things up a little.
+					*/
+      inline void set_active_cell(typename DH::active_cell_iterator &newcell);
   
-    /** Get ONE vector value at the given point. It is inefficient to
-	use single points. If you need more than one at a time, use the
-	vector_value_list function. For efficiency reasons, it is better
-	if all the points lie on the same cell. This is not mandatory,
-	however it does speed things up. */ 
-    virtual void vector_value (const Point<dim> &p,
-			       Vector<double>   &values) const;
+				       /**
+					* Get ONE vector value at the
+					* given point. It is
+					* inefficient to use single
+					* points. If you need more
+					* than one at a time, use the
+					* vector_value_list()
+					* function. For efficiency
+					* reasons, it is better if all
+					* the points lie on the same
+					* cell. This is not mandatory,
+					* however it does speed things
+					* up.
+					*/ 
+      virtual void vector_value (const Point<dim> &p,
+				 Vector<double>   &values) const;
 
-    /** Return the value of the function at the given point. Unless
-	there is only one component (i.e. the function is scalar), you
-	should state the component you want to have evaluated; it
-	defaults to zero, i.e. the first component.  It is inefficient
-	to use single points. If you need more than one at a time, use
-	the vector_value_list function. For efficiency reasons, it is
-	better if all the points lie on the same cell. This is not
-	mandatory, however it does speed things up. */
-    virtual double value (const Point< dim > &    p,
-			  const unsigned int  component = 0)    const;
+				       /**
+					* Return the value of the
+					* function at the given
+					* point. Unless there is only
+					* one component (i.e. the
+					* function is scalar), you
+					* should state the component
+					* you want to have evaluated;
+					* it defaults to zero,
+					* i.e. the first component.
+					* It is inefficient to use
+					* single points. If you need
+					* more than one at a time, use
+					* the vector_value_list
+					* function. For efficiency
+					* reasons, it is better if all
+					* the points lie on the same
+					* cell. This is not mandatory,
+					* however it does speed things
+					* up.
+					*/
+      virtual double value (const Point< dim > &    p,
+			    const unsigned int  component = 0)    const;
   
-    /** Set @p values to the point values of the specified component of
-	the function at the @p points. It is assumed that @p values
-	already has the right size, i.e. the same size as the points
-	array. This is rather efficient if all the points lie on the
-	same cell. If this is not the case, things may slow down a bit.
-    */
-    virtual void value_list (const std::vector<Point< dim > > &    points,
-			     std::vector< double > &values, 
-			     const unsigned int  component = 0)    const;
+				       /**
+					* Set @p values to the point
+					* values of the specified
+					* component of the function at
+					* the @p points. It is assumed
+					* that @p values already has
+					* the right size, i.e. the
+					* same size as the points
+					* array. This is rather
+					* efficient if all the points
+					* lie on the same cell. If
+					* this is not the case, things
+					* may slow down a bit.
+					*/
+      virtual void value_list (const std::vector<Point< dim > > &    points,
+			       std::vector< double > &values, 
+			       const unsigned int  component = 0)    const;
 
   
-    /** Set @p values to the point values of the function at the @p
-	points. It is assumed that @p values already has the right size,
-	i.e. the same size as the points array. This is rather efficient
-	if all the points lie on the same cell. If this is not the case,
-	things may slow down a bit.
-    */
-    virtual void vector_value_list (const std::vector<Point< dim > > &    points,
-				    std::vector< Vector<double> > &values) const;
+				       /**
+					* Set @p values to the point
+					* values of the function at
+					* the @p points. It is assumed
+					* that @p values already has
+					* the right size, i.e. the
+					* same size as the points
+					* array. This is rather
+					* efficient if all the points
+					* lie on the same cell. If
+					* this is not the case, things
+					* may slow down a bit.
+					*/
+      virtual void vector_value_list (const std::vector<Point< dim > > &    points,
+				      std::vector< Vector<double> > &values) const;
   
-    /** Return the gradient of all components of the function at the
-	given point.  It is inefficient to use single points. If you
-	need more than one at a time, use the vector_value_list
-	function. For efficiency reasons, it is better if all the points
-	lie on the same cell. This is not mandatory, however it does
-	speed things up. */
-    virtual void 
-    vector_gradient (const Point< dim > &p, 
-		     std::vector< Tensor< 1, dim > > &gradients) const;
+				       /**
+					* Return the gradient of all
+					* components of the function
+					* at the given point.  It is
+					* inefficient to use single
+					* points. If you need more
+					* than one at a time, use the
+					* vector_value_list
+					* function. For efficiency
+					* reasons, it is better if all
+					* the points lie on the same
+					* cell. This is not mandatory,
+					* however it does speed things
+					* up.
+					*/
+      virtual void 
+      vector_gradient (const Point< dim > &p, 
+		       std::vector< Tensor< 1, dim > > &gradients) const;
   
-    /** Return the gradient of the specified component of the function
-	at the given point.  It is inefficient to use single points. If
-	you need more than one at a time, use the vector_value_list
-	function. For efficiency reasons, it is better if all the points
-	lie on the same cell. This is not mandatory, however it does
-	speed things up. */
-    virtual Tensor<1,dim> gradient(const Point< dim > &p, 
-				   const unsigned int component = 0)const;
+				       /**
+					* Return the gradient of the
+					* specified component of the
+					* function at the given point.
+					* It is inefficient to use
+					* single points. If you need
+					* more than one at a time, use
+					* the vector_value_list
+					* function. For efficiency
+					* reasons, it is better if all
+					* the points lie on the same
+					* cell. This is not mandatory,
+					* however it does speed things
+					* up.
+					*/
+      virtual Tensor<1,dim> gradient(const Point< dim > &p, 
+				     const unsigned int component = 0)const;
   
-    /** Return the gradient of all components of the function at all
-	the given points.  This is rather efficient if all the points
-	lie on the same cell. If this is not the case, things may slow
-	down a bit. */
-    virtual void 
-    vector_gradient_list (const std::vector< Point< dim > > &p, 
-			  std::vector< 
-			  std::vector< Tensor< 1, dim > > > &gradients) const;
+				       /**
+					* Return the gradient of all
+					* components of the function
+					* at all the given points.
+					* This is rather efficient if
+					* all the points lie on the
+					* same cell. If this is not
+					* the case, things may slow
+					* down a bit.
+					*/
+      virtual void 
+      vector_gradient_list (const std::vector< Point< dim > > &p, 
+			    std::vector< 
+			    std::vector< Tensor< 1, dim > > > &gradients) const;
 
-    /** Return the gradient of the specified component of the function
-	at all the given points.  This is rather efficient if all the
-	points lie on the same cell. If this is not the case, things
-	may slow down a bit. */
-    virtual void 
-    gradient_list (const std::vector< Point< dim > > &p, 
-		   std::vector< Tensor< 1, dim > > &gradients, 
-		   const unsigned int component=0) const;
+				       /**
+					* Return the gradient of the
+					* specified component of the
+					* function at all the given
+					* points.  This is rather
+					* efficient if all the points
+					* lie on the same cell. If
+					* this is not the case, things
+					* may slow down a bit.
+					*/
+      virtual void 
+      gradient_list (const std::vector< Point< dim > > &p, 
+		     std::vector< Tensor< 1, dim > > &gradients, 
+		     const unsigned int component=0) const;
   
-    /** Create quadrature rules. This function groups the points into
-	blocks that live in the same cell, and fills up three vectors:
-	@p cells, @p qpoints and @p maps. The first is a list of the
-	cells that contain the points, the second is a list of
-	quadrature points matching each cell of the first list, and the
-	third contains the index of the given quadrature points, i.e.,
-	@p points[maps[3][4]] ends up as the 5th quadrature point in the
-	4th cell. This is where optimization would help. This function
-	returns the number of cells that contain the given set of
-	points.
-    */
-    unsigned int compute_point_locations(const std::vector<Point<dim> > &points,
-					 std::vector<typename DH::active_cell_iterator > &cells,
-					 std::vector<std::vector<Point<dim> > > &qpoints,
-					 std::vector<std::vector<unsigned int> > &maps) const;
+				       /**
+					* Create quadrature
+					* rules. This function groups
+					* the points into blocks that
+					* live in the same cell, and
+					* fills up three vectors: @p
+					* cells, @p qpoints and @p
+					* maps. The first is a list of
+					* the cells that contain the
+					* points, the second is a list
+					* of quadrature points
+					* matching each cell of the
+					* first list, and the third
+					* contains the index of the
+					* given quadrature points,
+					* i.e., @p points[maps[3][4]]
+					* ends up as the 5th
+					* quadrature point in the 4th
+					* cell. This is where
+					* optimization would
+					* help. This function returns
+					* the number of cells that
+					* contain the given set of
+					* points.
+					*/
+      unsigned int
+      compute_point_locations(const std::vector<Point<dim> > &points,
+			      std::vector<typename DH::active_cell_iterator > &cells,
+			      std::vector<std::vector<Point<dim> > > &qpoints,
+			      std::vector<std::vector<unsigned int> > &maps) const;
     
-  private:
-    /** Pointer to the dof handler. */
-    SmartPointer<const DH> dh;
+    private:
+				       /**
+					* Pointer to the dof handler.
+					*/
+      SmartPointer<const DH> dh;
     
-    /** A reference to the actual data vector. */
-    const VECTOR & data_vector;
+				       /**
+					* A reference to the actual
+					* data vector.
+					*/
+      const VECTOR & data_vector;
     
-    /** A reference to the mapping being used. */
-    const Mapping<dim> & mapping;
+				       /**
+					* A reference to the mapping
+					* being used.
+					*/
+      const Mapping<dim> & mapping;
   
-    /** The current cell in which we are evaluating*/
-    mutable typename DH::active_cell_iterator cell;
+				       /**
+					* The current cell in which we
+					* are evaluating.
+					*/
+      mutable typename DH::active_cell_iterator cell;
   
-    /** Store the number of components of this function. */
-    const unsigned int n_components;
+				       /**
+					* Store the number of
+					* components of this function.
+					*/
+      const unsigned int n_components;
   };
 }
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
diff --git a/deal.II/deal.II/include/numerics/fe_field_function.templates.h b/deal.II/deal.II/include/numerics/fe_field_function.templates.h
index 9f07fb7f55..9c66124947 100644
--- a/deal.II/deal.II/include/numerics/fe_field_function.templates.h
+++ b/deal.II/deal.II/include/numerics/fe_field_function.templates.h
@@ -11,34 +11,44 @@
 //
 //---------------------------------------------------------------------------
 
-#include <numerics/fe_field_function.h>
 #include <base/utilities.h>
 #include <base/logstream.h>
 #include <grid/grid_tools.h>
 #include <fe/fe_values.h>
+#include <numerics/fe_field_function.h>
+
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace Functions {
+namespace Functions
+{
 
   template <int dim, typename DH, typename VECTOR>
   FEFieldFunction<dim, DH, VECTOR>::FEFieldFunction (const DH &mydh, 
 						     const VECTOR &myv,
-						     const Mapping<dim> &mymapping) : 
-    Function<dim>(mydh.get_fe().n_components()),
-    dh(&mydh, "FEFieldFunction"),
-    data_vector(myv),
-    mapping(mymapping),
-    n_components(mydh.get_fe().n_components())
+						     const Mapping<dim> &mymapping)
+		  : 
+		  Function<dim>(mydh.get_fe().n_components()),
+		  dh(&mydh, "FEFieldFunction"),
+		  data_vector(myv),
+		  mapping(mymapping),
+		  n_components(mydh.get_fe().n_components())
   {
     cell = dh->begin_active();
   }
 
+
+  
   template <int dim, typename DH, typename VECTOR>
-  void FEFieldFunction<dim, DH, VECTOR>::set_active_cell(typename DH::active_cell_iterator &newcell) {
+  void
+  FEFieldFunction<dim, DH, VECTOR>::
+  set_active_cell(typename DH::active_cell_iterator &newcell)
+  {
     cell = newcell;
   }
-    
+
+
+  
   template <int dim, typename DH, typename VECTOR>
   void FEFieldFunction<dim, DH, VECTOR>::vector_value (const Point<dim> &p,
 						       Vector<double>   &values) const 
@@ -47,15 +57,16 @@ namespace Functions {
 	    ExcDimensionMismatch(values.size(), n_components));
     Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
   
-    // Check if we already have all we need
-    if(!GeometryInfo<dim>::is_inside_unit_cell(qp)) {
-      std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair 
-	= GridTools::find_active_cell_around_point (mapping, *dh, p); 
-      cell = my_pair.first;
-      qp = my_pair.second;
-    }				    
+				     // Check if we already have all we need
+    if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+      {
+	const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair 
+	  = GridTools::find_active_cell_around_point (mapping, *dh, p); 
+	cell = my_pair.first;
+	qp = my_pair.second;
+      }				    
   
-    // Now we can find out about the point
+				     // Now we can find out about the point
     Quadrature<dim> quad(qp);
     FEValues<dim> fe_v(mapping, dh->get_fe(), quad, 
 		       update_values);
@@ -65,9 +76,12 @@ namespace Functions {
     values = vvalues[0];
   }
 
+
+  
   template <int dim, typename DH, typename VECTOR>
-  double FEFieldFunction<dim, DH, VECTOR>::value 
-  (const Point<dim>   &p, unsigned int comp) const
+  double
+  FEFieldFunction<dim, DH, VECTOR>::value (const Point<dim>   &p,
+					   const unsigned int comp) const
   { 
     Vector<double> values(n_components);
     vector_value(p, values);
@@ -76,7 +90,8 @@ namespace Functions {
 
 
   template <int dim, typename DH, typename VECTOR>
-  void FEFieldFunction<dim, DH, VECTOR>::vector_gradient 
+  void
+  FEFieldFunction<dim, DH, VECTOR>::vector_gradient 
   (const Point<dim> &p,
    std::vector<Tensor<1,dim> > &gradients) const 
   { 
@@ -84,15 +99,16 @@ namespace Functions {
 	    ExcDimensionMismatch(gradients.size(), n_components));
     Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
   
-    // Check if we already have all we need
-    if(!GeometryInfo<dim>::is_inside_unit_cell(qp)) {
-      std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair  
-	= GridTools::find_active_cell_around_point (mapping, *dh, p); 
-      cell = my_pair.first;
-      qp = my_pair.second;
-    }				    
+				     // Check if we already have all we need
+    if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+      {
+	std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair  
+	  = GridTools::find_active_cell_around_point (mapping, *dh, p); 
+	cell = my_pair.first;
+	qp = my_pair.second;
+      }				    
   
-    // Now we can find out about the point
+				     // Now we can find out about the point
     Quadrature<dim> quad(qp);
     FEValues<dim> fe_v(mapping, dh->get_fe(), quad, 
 		       update_gradients);
@@ -103,6 +119,8 @@ namespace Functions {
     gradients = vgrads[0];
   }
 
+
+  
   template <int dim, typename DH, typename VECTOR>
   Tensor<1,dim> FEFieldFunction<dim, DH, VECTOR>::gradient 
   (const Point<dim>   &p, unsigned int comp) const
@@ -112,12 +130,14 @@ namespace Functions {
     return grads[comp];
   }
 
-  // Now the list versions
-  // ==============================
+				   // Now the list versions
+				   // ==============================
 
   template <int dim, typename DH, typename VECTOR>
-  void FEFieldFunction<dim, DH, VECTOR>::vector_value_list (const std::vector<Point< dim > > &    points,
-							    std::vector< Vector<double> > &values) const
+  void
+  FEFieldFunction<dim, DH, VECTOR>::
+  vector_value_list (const std::vector<Point< dim > > &    points,
+		     std::vector< Vector<double> > &values) const
   { 
     Assert(points.size() == values.size(),
 	   ExcDimensionMismatch(points.size(), values.size()));
@@ -128,44 +148,50 @@ namespace Functions {
   
     unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
   
-    // Now gather all the informations we need
-    for(unsigned int i=0; i<ncells; ++i) {
-      // Number of quadrature points on this cell
-      unsigned int nq = qpoints[i].size();
+				     // Now gather all the informations we need
+    for (unsigned int i=0; i<ncells; ++i)
+      {
+					 // Number of quadrature points on this cell
+	unsigned int nq = qpoints[i].size();
     
-      // Construct a quadrature formula
-      std::vector< double > ww(nq, 1./((double) nq));
-      Quadrature<dim> quad(qpoints[i], ww);
+					 // Construct a quadrature formula
+	std::vector< double > ww(nq, 1./((double) nq));
+	Quadrature<dim> quad(qpoints[i], ww);
     
-      // Get a function value object
-      FEValues<dim> fe_v(mapping, dh->get_fe(), quad, 
-			 update_values);
-      fe_v.reinit(cells[i]);
-      std::vector< Vector<double> > vvalues (nq, Vector<double>(n_components));
-      fe_v.get_function_values(data_vector, vvalues);
-      for(unsigned int q=0; q<nq; ++q)
-	values[maps[i][q]] = vvalues[q];
-    }
+					 // Get a function value object
+	FEValues<dim> fe_v(mapping, dh->get_fe(), quad, 
+			   update_values);
+	fe_v.reinit(cells[i]);
+	std::vector< Vector<double> > vvalues (nq, Vector<double>(n_components));
+	fe_v.get_function_values(data_vector, vvalues);
+	for (unsigned int q=0; q<nq; ++q)
+	  values[maps[i][q]] = vvalues[q];
+      }
   }
 
   template <int dim, typename DH, typename VECTOR>
-  void FEFieldFunction<dim, DH, VECTOR>::value_list (const std::vector<Point< dim > > &points,
-						     std::vector< double > &values, 
-						     const unsigned int  component) const
+  void
+  FEFieldFunction<dim, DH, VECTOR>::
+  value_list (const std::vector<Point< dim > > &points,
+	      std::vector< double > &values, 
+	      const unsigned int  component) const
   { 
     Assert(points.size() == values.size(),
 	   ExcDimensionMismatch(points.size(), values.size()));
     std::vector< Vector<double> > vvalues(points.size(), Vector<double>(n_components));
     vector_value_list(points, vvalues);
-    for(unsigned int q=0; q<points.size(); ++q)
+    for (unsigned int q=0; q<points.size(); ++q)
       values[q] = vvalues[q](component);
   }
 
 
+  
   template <int dim, typename DH, typename VECTOR>
-  void FEFieldFunction<dim, DH, VECTOR>::vector_gradient_list (const std::vector<Point< dim > > &    points,
-							       std::vector< 
-							       std::vector< Tensor<1,dim> > > &values) const
+  void
+  FEFieldFunction<dim, DH, VECTOR>::
+  vector_gradient_list (const std::vector<Point< dim > > &    points,
+			std::vector< 
+			std::vector< Tensor<1,dim> > > &values) const
   { 
     Assert(points.size() == values.size(),
 	   ExcDimensionMismatch(points.size(), values.size()));
@@ -176,131 +202,152 @@ namespace Functions {
   
     unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
   
-    // Now gather all the informations we need
-    for(unsigned int i=0; i<ncells; ++i) {
-      // Number of quadrature points on this cell
-      unsigned int nq = qpoints[i].size();
+				     // Now gather all the informations we need
+    for (unsigned int i=0; i<ncells; ++i)
+      {
+					 // Number of quadrature points on this cell
+	unsigned int nq = qpoints[i].size();
     
-      // Construct a quadrature formula
-      std::vector< double > ww(nq, 1./((double) nq));
-      Quadrature<dim> quad(qpoints[i], ww);
+					 // Construct a quadrature formula
+	std::vector< double > ww(nq, 1./((double) nq));
+	Quadrature<dim> quad(qpoints[i], ww);
     
-      // Get a function value object
-      FEValues<dim> fe_v(mapping, dh->get_fe(), quad, 
-			 update_gradients);
-      fe_v.reinit(cells[i]);
-      std::vector< std::vector<Tensor<1,dim> > > vgrads (nq, std::vector<Tensor<1,dim> >(n_components));
-      fe_v.get_function_grads(data_vector, vgrads);
-      for(unsigned int q=0; q<nq; ++q)
-	values[maps[i][q]] = vgrads[q];
-    }
+					 // Get a function value object
+	FEValues<dim> fe_v(mapping, dh->get_fe(), quad, 
+			   update_gradients);
+	fe_v.reinit(cells[i]);
+	std::vector< std::vector<Tensor<1,dim> > >
+	  vgrads (nq, std::vector<Tensor<1,dim> >(n_components));
+	fe_v.get_function_grads(data_vector, vgrads);
+	for (unsigned int q=0; q<nq; ++q)
+	  values[maps[i][q]] = vgrads[q];
+      }
   }
 
   template <int dim, typename DH, typename VECTOR>
-  void FEFieldFunction<dim, DH, VECTOR>::gradient_list (const std::vector<Point< dim > > &points,
-							std::vector< Tensor<1,dim> > &values, 
-							const unsigned int  component) const
+  void
+  FEFieldFunction<dim, DH, VECTOR>::
+  gradient_list (const std::vector<Point< dim > > &points,
+		 std::vector< Tensor<1,dim> > &values, 
+		 const unsigned int  component) const
   { 
     Assert(points.size() == values.size(),
 	   ExcDimensionMismatch(points.size(), values.size()));
-    std::vector< std::vector<Tensor<1,dim> > > vvalues(points.size(), std::vector<Tensor<1,dim> >(n_components));
+    std::vector< std::vector<Tensor<1,dim> > >
+      vvalues(points.size(), std::vector<Tensor<1,dim> >(n_components));
     vector_gradient_list(points, vvalues);
-    for(unsigned int q=0; q<points.size(); ++q)
+    for (unsigned int q=0; q<points.size(); ++q)
       values[q] = vvalues[q][component];
   }
 
-  // Now the distribute points function
+  
+
   template <int dim, typename DH, typename VECTOR>
   unsigned int FEFieldFunction<dim, DH, VECTOR>::
   compute_point_locations(const std::vector<Point<dim> > &points,
 			  std::vector<typename DH::active_cell_iterator > &cells,
 			  std::vector<std::vector<Point<dim> > > &qpoints,
-			  std::vector<std::vector<unsigned int> > &maps) const {
-    // How many points are here?
-    unsigned int np = points.size();
+			  std::vector<std::vector<unsigned int> > &maps) const
+  {
+				     // How many points are here?
+    const unsigned int np = points.size();
   
-    // Reset output maps.
+				     // Reset output maps.
     cells.clear();
     qpoints.clear();
     maps.clear();
     
-    // Now the easy case.
-    if(np==0) return 0;
+				     // Now the easy case.
+    if (np==0) return 0;
     
-    // Keep track of the points we found
+				     // Keep track of the points we
+				     // found
     std::vector<bool> point_flags(np, false);
   
-    // Set this to true untill all points have been classified
+				     // Set this to true untill all
+				     // points have been classified
     bool left_over = true;
   
-    // Current quadrature point
+				     // Current quadrature point
     Point<dim> qp = mapping.transform_real_to_unit_cell(cell, points[0]);
   
-    // Check if we already have a valid cell for the first point
-    if(!GeometryInfo<dim>::is_inside_unit_cell(qp)) {
-      std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair  = GridTools::find_active_cell_around_point
-	(mapping, *dh, points[0]); 
-      cell = my_pair.first;
-      qp = my_pair.second;
-      point_flags[0] = true;
-    }
+				     // Check if we already have a
+				     // valid cell for the first point
+    if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+      {
+	const std::pair<typename DH::active_cell_iterator, Point<dim> >
+	  my_pair  = GridTools::find_active_cell_around_point
+	  (mapping, *dh, points[0]); 
+	cell = my_pair.first;
+	qp = my_pair.second;
+	point_flags[0] = true;
+      }
     
-    // Put in the first point.
+				     // Put in the first point.
     cells.push_back(cell);
     qpoints.push_back(std::vector<Point<dim> >(1, qp));
     maps.push_back(std::vector<unsigned int> (1, 0));
   
-    // Check if we need to do anything else
-    if(points.size() > 1) {
+				     // Check if we need to do anything else
+    if (points.size() > 1)
       left_over = true;
-    } else {
+    else
       left_over = false;
-    }
+
   
-    // This is the first index of a non processed point
+				     // This is the first index of a non processed point
     unsigned int first_outside = 1;
   
-    // And this is the index of the current cell
+				     // And this is the index of the current cell
     unsigned int c = 0;
   
-    while(left_over == true) {
-      // Assume this is the last one
-      left_over = false;
-      Assert(first_outside < np,
-	     ExcIndexRange(first_outside, 0, np));
+    while (left_over == true)
+      {
+					 // Assume this is the last one
+	left_over = false;
+	Assert(first_outside < np,
+	       ExcIndexRange(first_outside, 0, np));
     
-      // If we found one in this cell, keep looking in the same cell
-      for(unsigned int p=first_outside; p<np; ++p) 
-	if(point_flags[p] == false) {
-	  Point<dim> qpoint =  mapping.transform_real_to_unit_cell(cell, points[p]);
-	  if(GeometryInfo<dim>::is_inside_unit_cell(qpoint)) {
-	    point_flags[p] = true;
-	    qpoints[c].push_back(qpoint);
-	    maps[c].push_back(p);
-	  } else {
-	    // Set things up for next round 
-	    if(left_over == false) first_outside = p;
-	    left_over = true;
+					 // If we found one in this cell, keep looking in the same cell
+	for (unsigned int p=first_outside; p<np; ++p) 
+	  if (point_flags[p] == false) {
+	    Point<dim> qpoint =  mapping.transform_real_to_unit_cell(cell, points[p]);
+	    if (GeometryInfo<dim>::is_inside_unit_cell(qpoint))
+	      {
+		point_flags[p] = true;
+		qpoints[c].push_back(qpoint);
+		maps[c].push_back(p);
+	      }
+	    else
+	      {
+						 // Set things up for next round 
+		if (left_over == false)
+		  first_outside = p;
+		left_over = true;
+	      }
 	  }
-	}
-      // If we got here and there is no left over, we are done. Else we
-      // need to find the next cell
-      if(left_over == true) {
-	std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair  
-	  = GridTools::find_active_cell_around_point (mapping, *dh, points[first_outside]); 
-	cells.push_back(my_pair.first);
-	qpoints.push_back(std::vector<Point<dim> >(1, my_pair.second));
-	maps.push_back(std::vector<unsigned int>(1, first_outside));
-	c++;
-	point_flags[first_outside] = true;
-	// And check if we can exit the loop now
-	if (first_outside == np-1) left_over = false;
-      }			
-    }
+					 // If we got here and there is
+					 // no left over, we are
+					 // done. Else we need to find
+					 // the next cell
+	if (left_over == true)
+	  {
+	    const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair  
+	      = GridTools::find_active_cell_around_point (mapping, *dh, points[first_outside]); 
+	    cells.push_back(my_pair.first);
+	    qpoints.push_back(std::vector<Point<dim> >(1, my_pair.second));
+	    maps.push_back(std::vector<unsigned int>(1, first_outside));
+	    c++;
+	    point_flags[first_outside] = true;
+					     // And check if we can exit the loop now
+	    if (first_outside == np-1)
+	      left_over = false;
+	  }			
+      }
   
-    // Augment of one the number of cells
+				     // Augment of one the number of cells
     ++c;
-    // Debug Checking
+				     // Debug Checking
     Assert(c == cells.size(), ExcInternalError());
   
     Assert(c == maps.size(),
@@ -311,17 +358,22 @@ namespace Functions {
   
 #ifdef DEBUG
     unsigned int qps = 0;
-    // The number of points in all the cells must be the same as the
-    // number of points we started off from.
-    for(unsigned int n=0; n<c; ++n) {
-      Assert(qpoints[n].size() == maps[n].size(),
-	     ExcDimensionMismatch(qpoints[n].size(), maps[n].size()));
-      qps += qpoints[n].size();
-    }
+				     // The number of points in all
+				     // the cells must be the same as
+				     // the number of points we
+				     // started off from.
+    for (unsigned int n=0; n<c; ++n)
+      {
+	Assert(qpoints[n].size() == maps[n].size(),
+	       ExcDimensionMismatch(qpoints[n].size(), maps[n].size()));
+	qps += qpoints[n].size();
+      }
     Assert(qps == np,
 	   ExcDimensionMismatch(qps, np));
 #endif
+    
     return c;
   }
 }
+
 DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/source/numerics/fe_field_function.cc b/deal.II/deal.II/source/numerics/fe_field_function.cc
index 9396dab1b7..d068d0b502 100644
--- a/deal.II/deal.II/source/numerics/fe_field_function.cc
+++ b/deal.II/deal.II/source/numerics/fe_field_function.cc
@@ -25,7 +25,8 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-namespace Functions {
+namespace Functions
+{
 
   template class FEFieldFunction<deal_II_dimension, 
 				 DoFHandler<deal_II_dimension>, 
@@ -40,8 +41,8 @@ namespace Functions {
 				 Vector<double> >;
 
   template class FEFieldFunction<deal_II_dimension, 
-			       MGDoFHandler<deal_II_dimension>, 
-			       BlockVector<double> >;
+				 MGDoFHandler<deal_II_dimension>, 
+				 BlockVector<double> >;
 
 #ifdef DEAL_II_USE_PETSC
 
@@ -50,27 +51,19 @@ namespace Functions {
 				 PETScWrappers::Vector >;
 
   template class FEFieldFunction<deal_II_dimension, 
-			       DoFHandler<deal_II_dimension>, 
-			       PETScWrappers::BlockVector >;
+				 DoFHandler<deal_II_dimension>, 
+				 PETScWrappers::BlockVector >;
 
   template class FEFieldFunction<deal_II_dimension, 
 				 MGDoFHandler<deal_II_dimension>, 
 				 PETScWrappers::Vector >;
 
   template class FEFieldFunction<deal_II_dimension, 
-			       MGDoFHandler<deal_II_dimension>, 
-			       PETScWrappers::BlockVector >;
+				 MGDoFHandler<deal_II_dimension>, 
+				 PETScWrappers::BlockVector >;
 
 #endif
   
-//   template class FEFieldFunction<deal_II_dimension, 
-// 				 DoFHandler<deal_II_dimension>, 
-// 				 Vector<double> >;
-
-//   template class FEFieldFunction<deal_II_dimension, 
-// 				 DoFHandler<deal_II_dimension>, 
-// 				 BlockVector<double> >;
-
 }
 
 DEAL_II_NAMESPACE_CLOSE
-- 
2.39.5