From 6e8b912c4ed34e33b0d7346bc6649423b8399e57 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 29 Aug 2007 21:46:46 +0000 Subject: [PATCH] Update comment and indentation style git-svn-id: https://svn.dealii.org/trunk@15092 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/parsed_function.h | 387 +++++++------ deal.II/base/source/parsed_function.cc | 113 ++-- .../include/numerics/fe_field_function.h | 508 +++++++++++------- .../numerics/fe_field_function.templates.h | 338 +++++++----- .../source/numerics/fe_field_function.cc | 23 +- 5 files changed, 807 insertions(+), 562 deletions(-) diff --git a/deal.II/base/include/base/parsed_function.h b/deal.II/base/include/base/parsed_function.h index 69da46e8dc..913a5db348 100644 --- a/deal.II/base/include/base/parsed_function.h +++ b/deal.II/base/include/base/parsed_function.h @@ -20,178 +20,241 @@ DEAL_II_NAMESPACE_OPEN -namespace Functions { - /** Friendly interface to the FunctionParser class. This class is - meant as a wrapper for the FunctionParser class. It provides two - methods to declare and parse a ParameterHandler object and creates - the Function object declared in the parameter file. This class is - derived from the AutoDerivativeFunction class, so you don't need - to specify derivatives. An example of usage of this class is as follows: - - \code - // A parameter handler - ParameterHandler prm; - - // Declare a section for the function we need - prm.enter_subsection("My vector function"); - ParsedFunction::declare_parameters(prm, dim); - prm.leave_subsection(); - - // Create a ParsedFunction - ParsedFunction 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::declare_parameters(prm, dim); + * prm.leave_subsection(); + * + * // Create a ParsedFunction + * ParsedFunction 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 class ParsedFunction : public AutoDerivativeFunction { - 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 &p, - Vector &values) const; + /** + * Get one value at the given + * point. + */ + virtual void vector_value (const Point &p, + Vector &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. */ - virtual void set_time(const double newtime); - private: - FunctionParser 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. + */ + virtual void set_time(const double newtime); + + private: + /** + * The object with which we do + * computations. + */ + FunctionParser function_object; }; } + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/base/source/parsed_function.cc b/deal.II/base/source/parsed_function.cc index 995eb64bb6..3ff828199c 100644 --- a/deal.II/base/source/parsed_function.cc +++ b/deal.II/base/source/parsed_function.cc @@ -17,45 +17,56 @@ DEAL_II_NAMESPACE_OPEN namespace Functions { template - ParsedFunction::ParsedFunction (const unsigned int n_components, const double h) : - AutoDerivativeFunction(h, n_components), - function_object(n_components) - { - } + ParsedFunction::ParsedFunction (const unsigned int n_components, const double h) + : + AutoDerivativeFunction(h, n_components), + function_object(n_components) + {} + + template - void ParsedFunction::declare_parameters(ParameterHandler &prm, unsigned int n_components) + void + ParsedFunction::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 void ParsedFunction::parse_parameters(ParameterHandler &prm) { @@ -66,33 +77,38 @@ namespace Functions { std::vector const_list = Utilities::split_string_list(constants_list, ','); std::map constants; - for(unsigned int i = 0; i < const_list.size(); ++i) { - std::vector 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 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 void ParsedFunction::vector_value (const Point &p, Vector &values) const @@ -100,6 +116,8 @@ namespace Functions { function_object.vector_value(p, values); } + + template double ParsedFunction::value (const Point &p, unsigned int comp) const @@ -107,6 +125,8 @@ namespace Functions { return function_object.value(p, comp); } + + template void ParsedFunction::set_time (const double newtime) { @@ -114,7 +134,8 @@ namespace Functions { AutoDerivativeFunction::set_time(newtime); } - // Explicit instantiations + +// Explicit instantiations template class ParsedFunction<1>; template class ParsedFunction<2>; template class ParsedFunction<3>; diff --git a/deal.II/deal.II/include/numerics/fe_field_function.h b/deal.II/deal.II/include/numerics/fe_field_function.h index 3137e18e31..61818e4fd3 100644 --- a/deal.II/deal.II/include/numerics/fe_field_function.h +++ b/deal.II/deal.II/include/numerics/fe_field_function.h @@ -25,226 +25,342 @@ DEAL_II_NAMESPACE_OPEN -namespace Functions { +namespace Functions +{ - /** This is an interpolation function for the given dof handler and - the given solution vector. The points at which this function can - be evaluated MUST be inside the domain of the dof handler, but - except from this, no other requirement is given. This function is - rather slow, as it needs to construct a quadrature object for the - point (or set of points) where you want to evaluate your finite - element function. In order to do so, it needs to find out where - the points lie. - - If you know in advance in which cell your points lye, you can - accelerate things a bit, by calling set_active_cell before - asking for values or gradients of the function. If you don't do - this, and your points don't lie in the cell that is currently - stored, the function GridTools::find_cell_around_point is called - to find out where the point is. You can specify an optional - mapping to use when looking for points in the grid. If you don't - do so, this function uses a Q1 mapping. - - Once the FEFieldFunction knows where the points lie, it creates a - quadrature formula for those points, and calls - FEValues::get_function_values or FEValues::get_function_grads with - the given quadrature points. - - If you only need the quadrature points but not the values of the - finite element function (you might want this for the adjoint - interpolation), you can also use the function @p - compute_point_locations alone. - - An example of how to use this function is the following: - - \code - - // Generate two triangulations - Triangulation tria_1; - Triangulation 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 dh1(tria_1); - Vector solution_1; - - // Do the same with the second - DoFHandler dh2; - Vector 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 fe_function_1 (dh_1, solution_1); - VectorTools::project(dh_2, constraints_2, quad, fe_function_1, solution_2); - - // Or interpolate it... - Vector 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 tria_1; + * Triangulation 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 dh1(tria_1); + * Vector solution_1; + * + * // Do the same with the second + * DoFHandler dh2; + * Vector 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 fe_function_1 (dh_1, solution_1); + * VectorTools::project(dh_2, constraints_2, quad, fe_function_1, solution_2); + * + * // Or interpolate it... + * Vector 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 , typename VECTOR=Vector > class FEFieldFunction : public Function { - 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 &mapping = StaticMappingQ1::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 &mapping = StaticMappingQ1::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 &p, - Vector &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 &p, + Vector &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 > & 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 > & 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 > & points, - std::vector< Vector > &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 > & points, + std::vector< Vector > &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 > &points, - std::vector &cells, - std::vector > > &qpoints, - std::vector > &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 > &points, + std::vector &cells, + std::vector > > &qpoints, + std::vector > &maps) const; - private: - /** Pointer to the dof handler. */ - SmartPointer dh; + private: + /** + * Pointer to the dof handler. + */ + SmartPointer 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 & mapping; + /** + * A reference to the mapping + * being used. + */ + const Mapping & mapping; - /** The current cell in which we are evaluating*/ - mutable typename DH::active_cell_iterator cell; + /** + * The current cell in which we + * are evaluating. + */ + mutable typename DH::active_cell_iterator cell; - /** Store the number of components of this function. */ - const unsigned int n_components; + /** + * Store the number of + * components of this function. + */ + const unsigned int n_components; }; } + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/deal.II/include/numerics/fe_field_function.templates.h b/deal.II/deal.II/include/numerics/fe_field_function.templates.h index 9f07fb7f55..9c66124947 100644 --- a/deal.II/deal.II/include/numerics/fe_field_function.templates.h +++ b/deal.II/deal.II/include/numerics/fe_field_function.templates.h @@ -11,34 +11,44 @@ // //--------------------------------------------------------------------------- -#include #include #include #include #include +#include + DEAL_II_NAMESPACE_OPEN -namespace Functions { +namespace Functions +{ template FEFieldFunction::FEFieldFunction (const DH &mydh, const VECTOR &myv, - const Mapping &mymapping) : - Function(mydh.get_fe().n_components()), - dh(&mydh, "FEFieldFunction"), - data_vector(myv), - mapping(mymapping), - n_components(mydh.get_fe().n_components()) + const Mapping &mymapping) + : + Function(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 - void FEFieldFunction::set_active_cell(typename DH::active_cell_iterator &newcell) { + void + FEFieldFunction:: + set_active_cell(typename DH::active_cell_iterator &newcell) + { cell = newcell; } - + + + template void FEFieldFunction::vector_value (const Point &p, Vector &values) const @@ -47,15 +57,16 @@ namespace Functions { ExcDimensionMismatch(values.size(), n_components)); Point qp = mapping.transform_real_to_unit_cell(cell, p); - // Check if we already have all we need - if(!GeometryInfo::is_inside_unit_cell(qp)) { - std::pair > 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::is_inside_unit_cell(qp)) + { + const std::pair > 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 quad(qp); FEValues fe_v(mapping, dh->get_fe(), quad, update_values); @@ -65,9 +76,12 @@ namespace Functions { values = vvalues[0]; } + + template - double FEFieldFunction::value - (const Point &p, unsigned int comp) const + double + FEFieldFunction::value (const Point &p, + const unsigned int comp) const { Vector values(n_components); vector_value(p, values); @@ -76,7 +90,8 @@ namespace Functions { template - void FEFieldFunction::vector_gradient + void + FEFieldFunction::vector_gradient (const Point &p, std::vector > &gradients) const { @@ -84,15 +99,16 @@ namespace Functions { ExcDimensionMismatch(gradients.size(), n_components)); Point qp = mapping.transform_real_to_unit_cell(cell, p); - // Check if we already have all we need - if(!GeometryInfo::is_inside_unit_cell(qp)) { - std::pair > 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::is_inside_unit_cell(qp)) + { + std::pair > 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 quad(qp); FEValues fe_v(mapping, dh->get_fe(), quad, update_gradients); @@ -103,6 +119,8 @@ namespace Functions { gradients = vgrads[0]; } + + template Tensor<1,dim> FEFieldFunction::gradient (const Point &p, unsigned int comp) const @@ -112,12 +130,14 @@ namespace Functions { return grads[comp]; } - // Now the list versions - // ============================== + // Now the list versions + // ============================== template - void FEFieldFunction::vector_value_list (const std::vector > & points, - std::vector< Vector > &values) const + void + FEFieldFunction:: + vector_value_list (const std::vector > & points, + std::vector< Vector > &values) const { Assert(points.size() == values.size(), ExcDimensionMismatch(points.size(), values.size())); @@ -128,44 +148,50 @@ namespace Functions { unsigned int ncells = compute_point_locations(points, cells, qpoints, maps); - // Now gather all the informations we need - for(unsigned int i=0; i ww(nq, 1./((double) nq)); - Quadrature quad(qpoints[i], ww); + // Construct a quadrature formula + std::vector< double > ww(nq, 1./((double) nq)); + Quadrature quad(qpoints[i], ww); - // Get a function value object - FEValues fe_v(mapping, dh->get_fe(), quad, - update_values); - fe_v.reinit(cells[i]); - std::vector< Vector > vvalues (nq, Vector(n_components)); - fe_v.get_function_values(data_vector, vvalues); - for(unsigned int q=0; q fe_v(mapping, dh->get_fe(), quad, + update_values); + fe_v.reinit(cells[i]); + std::vector< Vector > vvalues (nq, Vector(n_components)); + fe_v.get_function_values(data_vector, vvalues); + for (unsigned int q=0; q - void FEFieldFunction::value_list (const std::vector > &points, - std::vector< double > &values, - const unsigned int component) const + void + FEFieldFunction:: + value_list (const std::vector > &points, + std::vector< double > &values, + const unsigned int component) const { Assert(points.size() == values.size(), ExcDimensionMismatch(points.size(), values.size())); std::vector< Vector > vvalues(points.size(), Vector(n_components)); vector_value_list(points, vvalues); - for(unsigned int q=0; q - void FEFieldFunction::vector_gradient_list (const std::vector > & points, - std::vector< - std::vector< Tensor<1,dim> > > &values) const + void + FEFieldFunction:: + vector_gradient_list (const std::vector > & points, + std::vector< + std::vector< Tensor<1,dim> > > &values) const { Assert(points.size() == values.size(), ExcDimensionMismatch(points.size(), values.size())); @@ -176,131 +202,152 @@ namespace Functions { unsigned int ncells = compute_point_locations(points, cells, qpoints, maps); - // Now gather all the informations we need - for(unsigned int i=0; i ww(nq, 1./((double) nq)); - Quadrature quad(qpoints[i], ww); + // Construct a quadrature formula + std::vector< double > ww(nq, 1./((double) nq)); + Quadrature quad(qpoints[i], ww); - // Get a function value object - FEValues fe_v(mapping, dh->get_fe(), quad, - update_gradients); - fe_v.reinit(cells[i]); - std::vector< std::vector > > vgrads (nq, std::vector >(n_components)); - fe_v.get_function_grads(data_vector, vgrads); - for(unsigned int q=0; q fe_v(mapping, dh->get_fe(), quad, + update_gradients); + fe_v.reinit(cells[i]); + std::vector< std::vector > > + vgrads (nq, std::vector >(n_components)); + fe_v.get_function_grads(data_vector, vgrads); + for (unsigned int q=0; q - void FEFieldFunction::gradient_list (const std::vector > &points, - std::vector< Tensor<1,dim> > &values, - const unsigned int component) const + void + FEFieldFunction:: + gradient_list (const std::vector > &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 > > vvalues(points.size(), std::vector >(n_components)); + std::vector< std::vector > > + vvalues(points.size(), std::vector >(n_components)); vector_gradient_list(points, vvalues); - for(unsigned int q=0; q unsigned int FEFieldFunction:: compute_point_locations(const std::vector > &points, std::vector &cells, std::vector > > &qpoints, - std::vector > &maps) const { - // How many points are here? - unsigned int np = points.size(); + std::vector > &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 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 qp = mapping.transform_real_to_unit_cell(cell, points[0]); - // Check if we already have a valid cell for the first point - if(!GeometryInfo::is_inside_unit_cell(qp)) { - std::pair > 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::is_inside_unit_cell(qp)) + { + const std::pair > + 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 >(1, qp)); maps.push_back(std::vector (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 qpoint = mapping.transform_real_to_unit_cell(cell, points[p]); - if(GeometryInfo::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 qpoint = mapping.transform_real_to_unit_cell(cell, points[p]); + if (GeometryInfo::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 > my_pair - = GridTools::find_active_cell_around_point (mapping, *dh, points[first_outside]); - cells.push_back(my_pair.first); - qpoints.push_back(std::vector >(1, my_pair.second)); - maps.push_back(std::vector(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 > my_pair + = GridTools::find_active_cell_around_point (mapping, *dh, points[first_outside]); + cells.push_back(my_pair.first); + qpoints.push_back(std::vector >(1, my_pair.second)); + maps.push_back(std::vector(1, first_outside)); + c++; + point_flags[first_outside] = true; + // And check if we can exit the loop now + if (first_outside == np-1) + left_over = false; + } + } - // Augment of one the number of cells + // Augment of one the number of cells ++c; - // Debug Checking + // Debug Checking Assert(c == cells.size(), ExcInternalError()); Assert(c == maps.size(), @@ -311,17 +358,22 @@ namespace Functions { #ifdef DEBUG unsigned int qps = 0; - // The number of points in all the cells must be the same as the - // number of points we started off from. - for(unsigned int n=0; n, @@ -40,8 +41,8 @@ namespace Functions { Vector >; template class FEFieldFunction, - BlockVector >; + MGDoFHandler, + BlockVector >; #ifdef DEAL_II_USE_PETSC @@ -50,27 +51,19 @@ namespace Functions { PETScWrappers::Vector >; template class FEFieldFunction, - PETScWrappers::BlockVector >; + DoFHandler, + PETScWrappers::BlockVector >; template class FEFieldFunction, PETScWrappers::Vector >; template class FEFieldFunction, - PETScWrappers::BlockVector >; + MGDoFHandler, + PETScWrappers::BlockVector >; #endif -// template class FEFieldFunction, -// Vector >; - -// template class FEFieldFunction, -// BlockVector >; - } DEAL_II_NAMESPACE_CLOSE -- 2.39.5