--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
generalization of FilteredMatrix;
integration of the function parser library;
cubit journal file to export to ucd mesh format;
- FEFunction and ParsedFunction classes;
+ FEFieldFunction and ParsedFunction classes;
random bug fixes and enhancements.
<li><em>Oliver Kayser-Herold:</em>
</div>
</body>
-</html>
\ No newline at end of file
+</html>
<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.
<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
projects / dealii-svn.git / commitdiff