]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Update comment and indentation style
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 29 Aug 2007 21:46:46 +0000 (21:46 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 29 Aug 2007 21:46:46 +0000 (21:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@15092 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/parsed_function.h
deal.II/base/source/parsed_function.cc
deal.II/deal.II/include/numerics/fe_field_function.h
deal.II/deal.II/include/numerics/fe_field_function.templates.h
deal.II/deal.II/source/numerics/fe_field_function.cc

index 69da46e8dce4078effcf839e0d9e0601ff6e8e64..913a5db3488d19d4be7fde5e6301938a9dc0316e 100644 (file)
 
 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
index 995eb64bb6a676cd42b65db450c15c3dcc352981..3ff828199c6f2454d89beab49eced38b4dc402e0 100644 (file)
@@ -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>;
index 3137e18e310fd7ba8281a92eb5c8bf94dff3a415..61818e4fd31c79dca82b1395b399d9144b7edf9c 100644 (file)
 
 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
index 9f07fb7f557d0afd09f8a8a87a49bab3c43bdee6..9c66124947d47e1067d44bac0dcec39e058c66b4 100644 (file)
 //
 //---------------------------------------------------------------------------
 
-#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
index 9396dab1b7a4cc0236790de25de481abf14a6d65..d068d0b50270d1f77a78f4333ab1c0a1e85842f1 100644 (file)
@@ -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

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.