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
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)
{
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
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)
{
AutoDerivativeFunction<dim>::set_time(newtime);
}
- // Explicit instantiations
+
+// Explicit instantiations
template class ParsedFunction<1>;
template class ParsedFunction<2>;
template class ParsedFunction<3>;
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
//
//---------------------------------------------------------------------------
-#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
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);
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);
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
{
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);
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
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()));
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()));
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(),
#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
DEAL_II_NAMESPACE_OPEN
-namespace Functions {
+namespace Functions
+{
template class FEFieldFunction<deal_II_dimension,
DoFHandler<deal_II_dimension>,
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
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
projects / dealii-svn.git / commitdiff