]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added FEFieldFunction and ParsedFunction classes.
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 29 Aug 2007 20:16:16 +0000 (20:16 +0000)
committerLuca Heltai <luca.heltai@sissa.it>
Wed, 29 Aug 2007 20:16:16 +0000 (20:16 +0000)
git-svn-id: https://svn.dealii.org/trunk@15091 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/parsed_function.h [new file with mode: 0644]
deal.II/base/source/parsed_function.cc [new file with mode: 0644]
deal.II/deal.II/include/numerics/fe_field_function.h [new file with mode: 0644]
deal.II/deal.II/include/numerics/fe_field_function.templates.h [new file with mode: 0644]
deal.II/deal.II/source/numerics/fe_field_function.cc [new file with mode: 0644]
deal.II/doc/authors.html
deal.II/doc/news/changes.html

diff --git a/deal.II/base/include/base/parsed_function.h b/deal.II/base/include/base/parsed_function.h
new file mode 100644 (file)
index 0000000..69da46e
--- /dev/null
@@ -0,0 +1,197 @@
+//---------------------------------------------------------------------------
+//    $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __deal2__parsed_function_h
+#define __deal2__parsed_function_h
+
+#include <base/auto_derivative_function.h>
+#include <base/function_parser.h>
+#include <base/parameter_handler.h>
+
+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
+  */
+  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);
+  
+    /** 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);
+  
+    /** 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;
+  
+    /** 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;
+  };
+}
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/base/source/parsed_function.cc b/deal.II/base/source/parsed_function.cc
new file mode 100644 (file)
index 0000000..995eb64
--- /dev/null
@@ -0,0 +1,122 @@
+//---------------------------------------------------------------------------
+//    $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#include <base/parsed_function.h>
+#include <base/utilities.h>
+
+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)
+  {
+  }
+
+  template <int dim>
+  void ParsedFunction<dim>::declare_parameters(ParameterHandler &prm, 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;
+    }
+    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
+    std::string expr = "0";
+    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.");
+    prm.declare_entry("Function constants", "", Patterns::Anything(),
+                     "Any constant used inside the function which is not a variable name.");
+  }
+
+  template <int dim>
+  void ParsedFunction<dim>::parse_parameters(ParameterHandler &prm) 
+  {
+    std::string vnames = prm.get("Variable names");
+    std::string expression = prm.get("Function expression");
+    std::string constants_list = prm.get("Function constants");
+
+    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;
+    }
+    
+    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."));
+    }
+  }
+
+  template <int dim>
+  void ParsedFunction<dim>::vector_value (const Point<dim> &p,
+                                         Vector<double>   &values) const 
+  { 
+    function_object.vector_value(p, values);
+  }
+
+  template <int dim>
+  double ParsedFunction<dim>::value (const Point<dim>   &p,
+                                    unsigned int comp) const
+  { 
+    return function_object.value(p, comp);
+  }
+
+  template <int dim>
+  void ParsedFunction<dim>::set_time (const double newtime)
+  { 
+    function_object.set_time(newtime);
+    AutoDerivativeFunction<dim>::set_time(newtime);
+  }
+
+  // Explicit instantiations
+  template class ParsedFunction<1>;
+  template class ParsedFunction<2>;
+  template class ParsedFunction<3>;
+}
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/include/numerics/fe_field_function.h b/deal.II/deal.II/include/numerics/fe_field_function.h
new file mode 100644 (file)
index 0000000..3137e18
--- /dev/null
@@ -0,0 +1,250 @@
+//---------------------------------------------------------------------------
+//    $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__fe_function_h
+#define __deal2__fe_function_h
+
+#include <base/function.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/mapping_q1.h>
+#include <base/function.h>
+#include <base/point.h>
+#include <base/tensor.h>
+
+#include <lac/vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+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
+  */
+  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);
+  
+    /** 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;
+
+    /** 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 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 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 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;
+    
+  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 mapping being used. */
+    const Mapping<dim> & mapping;
+  
+    /** 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;
+  };
+}
+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
new file mode 100644 (file)
index 0000000..9f07fb7
--- /dev/null
@@ -0,0 +1,327 @@
+//---------------------------------------------------------------------------
+//    $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/fe_field_function.h>
+#include <base/utilities.h>
+#include <base/logstream.h>
+#include <grid/grid_tools.h>
+#include <fe/fe_values.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+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())
+  {
+    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) {
+    cell = newcell;
+  }
+    
+  template <int dim, typename DH, typename VECTOR>
+  void FEFieldFunction<dim, DH, VECTOR>::vector_value (const Point<dim> &p,
+                                                      Vector<double>   &values) const 
+  { 
+    Assert (values.size() == n_components, 
+           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;
+    }                              
+  
+    // Now we can find out about the point
+    Quadrature<dim> quad(qp);
+    FEValues<dim> fe_v(mapping, dh->get_fe(), quad, 
+                      update_values);
+    fe_v.reinit(cell);
+    std::vector< Vector<double> > vvalues (1, values);
+    fe_v.get_function_values(data_vector, vvalues);
+    values = vvalues[0];
+  }
+
+  template <int dim, typename DH, typename VECTOR>
+  double FEFieldFunction<dim, DH, VECTOR>::value 
+  (const Point<dim>   &p, unsigned int comp) const
+  { 
+    Vector<double> values(n_components);
+    vector_value(p, values);
+    return values(comp);
+  }
+
+
+  template <int dim, typename DH, typename VECTOR>
+  void FEFieldFunction<dim, DH, VECTOR>::vector_gradient 
+  (const Point<dim> &p,
+   std::vector<Tensor<1,dim> > &gradients) const 
+  { 
+    Assert (gradients.size() == n_components, 
+           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;
+    }                              
+  
+    // Now we can find out about the point
+    Quadrature<dim> quad(qp);
+    FEValues<dim> fe_v(mapping, dh->get_fe(), quad, 
+                      update_gradients);
+    fe_v.reinit(cell);
+    std::vector< std::vector<Tensor<1,dim> > > vgrads 
+      (1,  std::vector<Tensor<1,dim> >(n_components) );
+    fe_v.get_function_grads(data_vector, vgrads);
+    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
+  { 
+    std::vector<Tensor<1,dim> > grads(n_components);
+    vector_gradient(p, grads);
+    return grads[comp];
+  }
+
+  // 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
+  { 
+    Assert(points.size() == values.size(),
+          ExcDimensionMismatch(points.size(), values.size()));
+
+    std::vector<typename DH::active_cell_iterator > cells;
+    std::vector<std::vector<Point<dim> > > qpoints;
+    std::vector<std::vector<unsigned int> > maps;
+  
+    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();
+    
+      // 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];
+    }
+  }
+
+  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
+  { 
+    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)
+      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
+  { 
+    Assert(points.size() == values.size(),
+          ExcDimensionMismatch(points.size(), values.size()));
+
+    std::vector<typename DH::active_cell_iterator > cells;
+    std::vector<std::vector<Point<dim> > > qpoints;
+    std::vector<std::vector<unsigned int> > maps;
+  
+    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();
+    
+      // 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];
+    }
+  }
+
+  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
+  { 
+    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));
+    vector_gradient_list(points, vvalues);
+    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();
+  
+    // Reset output maps.
+    cells.clear();
+    qpoints.clear();
+    maps.clear();
+    
+    // Now the easy case.
+    if(np==0) return 0;
+    
+    // Keep track of the points we found
+    std::vector<bool> point_flags(np, false);
+  
+    // Set this to true untill all points have been classified
+    bool left_over = true;
+  
+    // 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;
+    }
+    
+    // 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) {
+      left_over = true;
+    } else {
+      left_over = false;
+    }
+  
+    // This is the first index of a non processed point
+    unsigned int first_outside = 1;
+  
+    // 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));
+    
+      // 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;
+      }                        
+    }
+  
+    // Augment of one the number of cells
+    ++c;
+    // Debug Checking
+    Assert(c == cells.size(), ExcInternalError());
+  
+    Assert(c == maps.size(),
+          ExcDimensionMismatch(c, maps.size()));
+  
+    Assert(c == qpoints.size(),
+          ExcDimensionMismatch(c, qpoints.size()));
+  
+#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();
+    }
+    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
new file mode 100644 (file)
index 0000000..9396dab
--- /dev/null
@@ -0,0 +1,76 @@
+//---------------------------------------------------------------------------
+//    $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/fe_field_function.templates.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <hp/dof_handler.h>
+
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/petsc_vector.h>
+#include <lac/petsc_block_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Functions {
+
+  template class FEFieldFunction<deal_II_dimension, 
+                                DoFHandler<deal_II_dimension>, 
+                                Vector<double> >;
+
+  template class FEFieldFunction<deal_II_dimension, 
+                                DoFHandler<deal_II_dimension>, 
+                                BlockVector<double> >;
+
+  template class FEFieldFunction<deal_II_dimension, 
+                                MGDoFHandler<deal_II_dimension>, 
+                                Vector<double> >;
+
+  template class FEFieldFunction<deal_II_dimension, 
+                              MGDoFHandler<deal_II_dimension>, 
+                              BlockVector<double> >;
+
+#ifdef DEAL_II_USE_PETSC
+
+  template class FEFieldFunction<deal_II_dimension, 
+                                DoFHandler<deal_II_dimension>, 
+                                PETScWrappers::Vector >;
+
+  template class FEFieldFunction<deal_II_dimension, 
+                              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 >;
+
+#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
index e4a40a724a63f3666a36822d9f2a4c53caa3a93c..cf24a1e525513b196a13a81b290f29211078e216 100644 (file)
@@ -69,7 +69,7 @@ alphabetical order):
     generalization of FilteredMatrix; &nbsp;
     integration of the function parser library; &nbsp;
     cubit journal file to export to ucd mesh format; &nbsp;
-    FEFunction and ParsedFunction classes; &nbsp;
+    FEFieldFunction and ParsedFunction classes; &nbsp;
     random bug fixes and enhancements.
 
 <li><em>Oliver Kayser-Herold:</em>
@@ -157,4 +157,4 @@ alphabetical order):
  </div>
 
 </body>
-</html>
\ No newline at end of file
+</html>
index 7c15ac8058a074ad2b08583081f71a3051f8d9fe..8de754400803b9f7e8433169bf28a79d1dbf941a 100644 (file)
@@ -481,6 +481,13 @@ inconvenience this causes.
 
 <ol>
 
+  <li> <p> New: There is a new class:
+       <code>Functions::FEFieldFunction</code> which is a Function
+       interface to a finite element solution. 
+       <br> 
+       (Luca Heltai 2007/08/29)
+
+  
   <li> <p> Improved: <code class="class">FunctionDerivative</code> is now
   derived from <code class="class">AutoDerivativeFunction</code> and implements
   gradients as well, giving you automatic second derivatives of a function.
@@ -1115,6 +1122,12 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p> New: There is a new class:
+       <code>Functions::ParsedFunction</code> which is friendly
+       wrapper to the <code>FunctionParser</code> class.
+       <br> 
+       (Luca Heltai 2007/08/29)
+  
    <li> <p>Fixed: the function
        <code>DataOut::build_patches</code> 
        had a quadratic algorithm when generatic cell-data (as opposed

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.