From 5b17d4c01f638b5886d44433cb7c89c0ebcbcf6e Mon Sep 17 00:00:00 2001 From: Konrad Simon Date: Mon, 16 Sep 2019 23:15:01 +0200 Subject: [PATCH] New feature class TensorFunctionParser. The class inherits directly from TensorFuntion and has an interface similar to FunctionParser. Some stuff to iterface muParser was moved to different files mu_parser_internal.h and *.cc. The documentation provided is self contained apart from some places where we refer to the documentation of FuntionParser. Changelog file provided. (All commits are squashed into this one.) --- doc/news/changes/minor/20191217KonradSimon | 4 + include/deal.II/base/mu_parser_internal.h | 97 +++++ include/deal.II/base/tensor_function_parser.h | 361 ++++++++++++++++ source/base/CMakeLists.txt | 3 + source/base/function_parser.cc | 225 ++-------- source/base/mu_parser_internal.cc | 203 +++++++++ source/base/tensor_function_parser.cc | 400 ++++++++++++++++++ source/base/tensor_function_parser.inst.in | 25 ++ tests/base/tensor_function_parser_01.cc | 69 +++ ...nction_parser_01.with_muparser=true.output | 2 + tests/base/tensor_function_parser_02.cc | 69 +++ ...nction_parser_02.with_muparser=true.output | 2 + tests/base/tensor_function_parser_03.cc | 82 ++++ ...nction_parser_03.with_muparser=true.output | 17 + tests/base/tensor_function_parser_04.cc | 78 ++++ ...nction_parser_04.with_muparser=true.output | 12 + 16 files changed, 1468 insertions(+), 181 deletions(-) create mode 100644 doc/news/changes/minor/20191217KonradSimon create mode 100644 include/deal.II/base/mu_parser_internal.h create mode 100644 include/deal.II/base/tensor_function_parser.h create mode 100644 source/base/mu_parser_internal.cc create mode 100644 source/base/tensor_function_parser.cc create mode 100644 source/base/tensor_function_parser.inst.in create mode 100644 tests/base/tensor_function_parser_01.cc create mode 100644 tests/base/tensor_function_parser_01.with_muparser=true.output create mode 100644 tests/base/tensor_function_parser_02.cc create mode 100644 tests/base/tensor_function_parser_02.with_muparser=true.output create mode 100644 tests/base/tensor_function_parser_03.cc create mode 100644 tests/base/tensor_function_parser_03.with_muparser=true.output create mode 100644 tests/base/tensor_function_parser_04.cc create mode 100644 tests/base/tensor_function_parser_04.with_muparser=true.output diff --git a/doc/news/changes/minor/20191217KonradSimon b/doc/news/changes/minor/20191217KonradSimon new file mode 100644 index 0000000000..05f2c048cc --- /dev/null +++ b/doc/news/changes/minor/20191217KonradSimon @@ -0,0 +1,4 @@ +New: The TensorFunctionParser class allows to read +in a function similarly to the FunctionParser class using the MuParser. +
+(Konrad Simon, 2019/12/17) \ No newline at end of file diff --git a/include/deal.II/base/mu_parser_internal.h b/include/deal.II/base/mu_parser_internal.h new file mode 100644 index 0000000000..00e6d95f01 --- /dev/null +++ b/include/deal.II/base/mu_parser_internal.h @@ -0,0 +1,97 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_mu_parser_internal_h +#define dealii_mu_parser_internal_h + +// This file contains functions used internally by the FunctionParser +// and the TensorFunctionParser class. + +#include + +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + + +#ifdef DEAL_II_WITH_MUPARSER + +namespace internal +{ + namespace FunctionParser + { + int + mu_round(double val); + + double + mu_if(double condition, double thenvalue, double elsevalue); + + double + mu_or(double left, double right); + + double + mu_and(double left, double right); + + double + mu_int(double value); + + double + mu_ceil(double value); + + double + mu_floor(double value); + + double + mu_cot(double value); + + double + mu_csc(double value); + + double + mu_sec(double value); + + double + mu_log(double value); + + double + mu_pow(double a, double b); + + double + mu_erfc(double value); + + // returns a random value in the range [0,1] initializing the generator + // with the given seed + double + mu_rand_seed(double seed); + + // returns a random value in the range [0,1] + double + mu_rand(); + + extern std::vector function_names; + + } // namespace FunctionParser + +} // namespace internal +#endif + + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/base/tensor_function_parser.h b/include/deal.II/base/tensor_function_parser.h new file mode 100644 index 0000000000..4137b5198b --- /dev/null +++ b/include/deal.II/base/tensor_function_parser.h @@ -0,0 +1,361 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_tensor_function_parser_h +#define dealii_tensor_function_parser_h + + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace mu +{ + class Parser; +} + +DEAL_II_NAMESPACE_OPEN + +// Forward declaration +#ifndef DOXYGEN +template +class Vector; +#endif + +/** + * This class implements a tensor function object that gets its value by parsing + * a string describing this function. It is a wrapper class for the muparser + * library (see http://muparser.beltoforion.de/). This class is essentially an + * extension of the FunctionParser class to read in a TensorFunction. The class + * reads in an expression of length dimrank (separated by a + * semicolon) where the components of the tensor function are filled according + * to the C++ convention (fastest index is the most right one). + * + * @note In contrast to the FunctionParser class the TensorFunctionParser class does not support + * automatic differentiation. + * + * A minimal example for the usage of the class would be: + * @code + * // set up time dependent tensor function: + * const std::string variables = "x,y,t"; + * const std::string expression = + * "exp(-t)*cos(x+y);-sin(pi*x*y-t);sin(pi*x*y-t);exp(t)*cos(x+y)"; + * std::map constants; + * constants["pi"] = numbers::PI; + * + * // TensorFunctionParser with 2+1 variables (space + time) in 2D of rank 2. + * // It is necessary to tell the parser that there is an additional variable + * // to be taken into account (t). + * TensorFunctionParser<2,2> tfp; + * tfp.initialize(variables, + * expression, + * constants, + * true); // flag for time dependence + * + * // Point at which we want to evaluate the function + * Point<2> point(0.0, 1.0); + * + * // evaluate the expression at 'point': + * double result = tfp.value(point); + * + * deallog << "Function '" << expression << "'" + * << " @ " << point + * << " is: " + * << std::endl + * << result[0][0] << " " << result[0][1] << std::endl + * << result[1][0] << " " << result[1][1] + * << std::endl; + * @endcode + * + * See also the documentation of the FunctionParser class. + * + * This class overloads the virtual method value() and value_list() of the + * TensorFunction base class with the byte compiled versions of the expressions + * given to the initialize() methods. Note that the class will not work unless + * you first call the initialize() method that accepts the text description of + * the function as an argument (among other things). + * + * The syntax to describe a function follows usual programming practice, and + * is explained in detail at the homepage of the underlying muparser library + * at http://muparser.beltoforion.de/ . + * + * + * Vector-valued functions can either be declared using strings where the + * function components are separated by semicolons, or using a vector of + * strings each defining one vector component. + * + * + * @ingroup functions + * @author Konrad Simon, 2019 + */ +template +class TensorFunctionParser : public TensorFunction +{ +public: + /** + * Standard constructor. Only set initial time. This object needs to be + * initialized with the initialize() method before you can use it. If an + * attempt to use this function is made before the initialize() method has + * been called, then an exception is thrown. + */ + TensorFunctionParser(const double initial_time = 0.0); + + /** + * Constructor for parsed functions. This object needs to be initialized + * with the initialize() method before you can use it. If an attempt to + * use this function is made before the initialize() method has been called, + * then an exception is thrown. + * Takes a semicolon separated list of expressions (one for each component + * of the tensor function), an optional comma-separated list of constants. + */ + TensorFunctionParser( + const std::string &expression, + const std::string &constants = "", + const std::string &variable_names = default_variable_names() + ",t"); + + /** + * Copy constructor. Objects of this type can not be copied, and + * consequently this constructor is deleted. + */ + TensorFunctionParser(const TensorFunctionParser &) = delete; + + /** + * Move constructor. Objects of this type can not be moved, and + * consequently this constructor is deleted. + */ + TensorFunctionParser(TensorFunctionParser &&) = delete; + + /** + * Destructor. + */ + virtual ~TensorFunctionParser() override; + + /** + * Copy operator. Objects of this type can not be copied, and + * consequently this operator is deleted. + */ + TensorFunctionParser & + operator=(const TensorFunctionParser &) = delete; + + /** + * Move operator. Objects of this type can not be moved, and + * consequently this operator is deleted. + */ + TensorFunctionParser & + operator=(TensorFunctionParser &&) = delete; + + /** + * Type for the constant map. Used by the initialize() method. + */ + using ConstMap = std::map; + + + /** + * Initialize the tensor function. This method accepts the following + * parameters: + * + * @param[in] vars A string with the variables that will be used by the + * expressions to be evaluated. Note that the variables can have any name + * (of course different from the function names defined above!), but the + * order IS important. The first variable will correspond to the first + * component of the point in which the function is evaluated, the second + * variable to the second component and so forth. If this function is also + * time dependent, then it is necessary to specify it by setting the + * time_dependent parameter to true. An exception is thrown if + * the number of variables specified here is different from dim (if this + * function is not time-dependent) or from dim+1 (if it is time-dependent). + * + * @param[in] expressions A vector of strings containing the expressions that + * will be byte compiled by the internal parser (TensorFunctionParser). Note + * that the size of this vector must match exactly the number of components of + * the TensorFunctionParser, as declared in the constructor. If this is not + * the case, an exception is thrown. + * + * @param[in] constants A map of constants used to pass any necessary constant + * that we want to specify in our expressions (in the example above the + * number pi). An expression is valid if and only if it contains only + * defined variables and defined constants (other than the functions + * specified above). If a constant is given whose name is not valid (eg: + * constants["sin"] = 1.5;) an exception is thrown. + * + * @param[in] time_dependent If this is a time dependent function, then the + * last variable declared in vars is assumed to be the time variable, + * and this->get_time() is used to initialize it when evaluating the + * function. Naturally the number of variables parsed by the initialize() + * method in this case is dim+1. The value of this parameter defaults to + * false, i.e. do not consider time. + */ + void + initialize(const std::string & vars, + const std::vector &expressions, + const ConstMap & constants, + const bool time_dependent = false); + + /** + * Initialize the function. Same as above, but accepts a string rather than + * a vector of strings. If this is a vector valued function, its components + * are expected to be separated by a semicolon. An exception is thrown if + * this method is called and the number of components successfully parsed + * does not match the number of components of the base function. + */ + void + initialize(const std::string &vars, + const std::string &expression, + const ConstMap & constants, + const bool time_dependent = false); + + /** + * A function that returns default names for variables, to be used in the + * first argument of the initialize() functions: it returns "x" in 1d, "x,y" + * in 2d, and "x,y,z" in 3d. + */ + static std::string + default_variable_names(); + + /** + * Return the value of the tensor function at the given point. + */ + virtual Tensor + value(const Point &p) const override; + + /** + * Return the value of the tensor function at the given point. + */ + virtual void + value_list(const std::vector> & p, + std::vector> &values) const override; + + /** + * Return an array of function expressions (one per component), used to + * initialize this function. + */ + const std::vector & + get_expressions() const; + + /** + * @addtogroup Exceptions + * @{ + */ + DeclException2(ExcParseError, + int, + std::string, + << "Parsing Error at Column " << arg1 + << ". The parser said: " << arg2); + + DeclException2(ExcInvalidExpressionSize, + int, + int, + << "The number of components (" << arg1 + << ") is not equal to the number of expressions (" << arg2 + << ")."); + + //@} + +private: +#ifdef DEAL_II_WITH_MUPARSER + /** + * Place for the variables for each thread + */ + mutable Threads::ThreadLocalStorage> vars; + + /** + * The muParser objects for each thread (and one for each component). We are + * storing a unique_ptr so that we don't need to include the definition of + * mu::Parser in this header. + */ + mutable Threads::ThreadLocalStorage>> + tfp; + + /** + * An array to keep track of all the constants, required to initialize tfp in + * each thread. + */ + std::map constants; + + /** + * An array for the variable names, required to initialize tfp in each + * thread. + */ + std::vector var_names; + + /** + * Initialize tfp and vars on the current thread. This function may only be + * called once per thread. A thread can test whether the function has + * already been called by testing whether 'tfp.get().size()==0' (not + * initialized) or >0 (already initialized). + */ + void + init_muparser() const; +#endif + + /** + * An array of function expressions (one per component), required to + * initialize tfp in each thread. + */ + std::vector expressions; + + /** + * State of usability. This variable is checked every time the function is + * called for evaluation. It's set to true in the initialize() methods. + */ + bool initialized; + + /** + * Number of variables. If this is also a function of time, then the number + * of variables is dim+1, otherwise it is dim. In the case that this is a + * time dependent function, the time is supposed to be the last variable. If + * #n_vars is not identical to the number of the variables parsed by the + * initialize() method, then an exception is thrown. + */ + unsigned int n_vars; + + /** + * Number of components is equal dimrank<\sup>. + */ + unsigned int n_components; +}; + + +template +std::string +TensorFunctionParser::default_variable_names() +{ + switch (dim) + { + case 1: + return "x"; + case 2: + return "x,y"; + case 3: + return "x,y,z"; + default: + Assert(false, ExcNotImplemented()); + } + return ""; +} + + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 6920744d5f..fbf5754950 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -45,6 +45,7 @@ SET(_unity_include_src logstream.cc hdf5.cc mpi.cc + mu_parser_internal.cc multithread_info.cc named_selection.cc numbers.cc @@ -79,6 +80,7 @@ SET(_unity_include_src subscriptor.cc table_handler.cc tensor_function.cc + tensor_function_parser.cc tensor_product_polynomials.cc tensor_product_polynomials_bubbles.cc tensor_product_polynomials_const.cc @@ -128,6 +130,7 @@ SET(_inst symmetric_tensor.inst.in tensor.inst.in tensor_function.inst.in + tensor_function_parser.inst.in time_stepping.inst.in ) diff --git a/source/base/function_parser.cc b/source/base/function_parser.cc index 1a0b4d33bf..0a20249014 100644 --- a/source/base/function_parser.cc +++ b/source/base/function_parser.cc @@ -15,6 +15,7 @@ #include +#include #include #include #include @@ -133,130 +134,6 @@ FunctionParser::initialize(const std::string & variables, -namespace internal -{ - // convert double into int - int - mu_round(double val) - { - return static_cast(val + ((val >= 0.0) ? 0.5 : -0.5)); - } - - double - mu_if(double condition, double thenvalue, double elsevalue) - { - if (mu_round(condition)) - return thenvalue; - else - return elsevalue; - } - - double - mu_or(double left, double right) - { - return (mu_round(left)) || (mu_round(right)); - } - - double - mu_and(double left, double right) - { - return (mu_round(left)) && (mu_round(right)); - } - - double - mu_int(double value) - { - return static_cast(mu_round(value)); - } - - double - mu_ceil(double value) - { - return std::ceil(value); - } - - double - mu_floor(double value) - { - return std::floor(value); - } - - double - mu_cot(double value) - { - return 1.0 / std::tan(value); - } - - double - mu_csc(double value) - { - return 1.0 / std::sin(value); - } - - double - mu_sec(double value) - { - return 1.0 / std::cos(value); - } - - double - mu_log(double value) - { - return std::log(value); - } - - double - mu_pow(double a, double b) - { - return std::pow(a, b); - } - - double - mu_erfc(double value) - { - return std::erfc(value); - } - - // returns a random value in the range [0,1] initializing the generator - // with the given seed - double - mu_rand_seed(double seed) - { - static Threads::Mutex rand_mutex; - std::lock_guard lock(rand_mutex); - - static boost::random::uniform_real_distribution<> uniform_distribution(0, - 1); - - // for each seed an unique random number generator is created, - // which is initialized with the seed itself - // we could use std::mt19937 but doing so results in compiler-dependent - // output. - static std::map rng_map; - - if (rng_map.find(seed) == rng_map.end()) - rng_map[seed] = boost::random::mt19937(static_cast(seed)); - - return uniform_distribution(rng_map[seed]); - } - - // returns a random value in the range [0,1] - double - mu_rand() - { - static Threads::Mutex rand_mutex; - std::lock_guard lock(rand_mutex); - static boost::random::uniform_real_distribution<> uniform_distribution(0, - 1); - static boost::random::mt19937 rng( - static_cast(std::time(nullptr))); - return uniform_distribution(rng); - } - -} // namespace internal - - - template void FunctionParser::init_muparser() const @@ -283,20 +160,44 @@ FunctionParser::init_muparser() const fp.get()[component]->DefineVar(var_names[iv], &vars.get()[iv]); // define some compatibility functions: - fp.get()[component]->DefineFun("if", internal::mu_if, true); - fp.get()[component]->DefineOprt("|", internal::mu_or, 1); - fp.get()[component]->DefineOprt("&", internal::mu_and, 2); - fp.get()[component]->DefineFun("int", internal::mu_int, true); - fp.get()[component]->DefineFun("ceil", internal::mu_ceil, true); - fp.get()[component]->DefineFun("cot", internal::mu_cot, true); - fp.get()[component]->DefineFun("csc", internal::mu_csc, true); - fp.get()[component]->DefineFun("floor", internal::mu_floor, true); - fp.get()[component]->DefineFun("sec", internal::mu_sec, true); - fp.get()[component]->DefineFun("log", internal::mu_log, true); - fp.get()[component]->DefineFun("pow", internal::mu_pow, true); - fp.get()[component]->DefineFun("erfc", internal::mu_erfc, true); - fp.get()[component]->DefineFun("rand_seed", internal::mu_rand_seed, true); - fp.get()[component]->DefineFun("rand", internal::mu_rand, true); + fp.get()[component]->DefineFun("if", + internal::FunctionParser::mu_if, + true); + fp.get()[component]->DefineOprt("|", internal::FunctionParser::mu_or, 1); + fp.get()[component]->DefineOprt("&", internal::FunctionParser::mu_and, 2); + fp.get()[component]->DefineFun("int", + internal::FunctionParser::mu_int, + true); + fp.get()[component]->DefineFun("ceil", + internal::FunctionParser::mu_ceil, + true); + fp.get()[component]->DefineFun("cot", + internal::FunctionParser::mu_cot, + true); + fp.get()[component]->DefineFun("csc", + internal::FunctionParser::mu_csc, + true); + fp.get()[component]->DefineFun("floor", + internal::FunctionParser::mu_floor, + true); + fp.get()[component]->DefineFun("sec", + internal::FunctionParser::mu_sec, + true); + fp.get()[component]->DefineFun("log", + internal::FunctionParser::mu_log, + true); + fp.get()[component]->DefineFun("pow", + internal::FunctionParser::mu_pow, + true); + fp.get()[component]->DefineFun("erfc", + internal::FunctionParser::mu_erfc, + true); + fp.get()[component]->DefineFun("rand_seed", + internal::FunctionParser::mu_rand_seed, + true); + fp.get()[component]->DefineFun("rand", + internal::FunctionParser::mu_rand, + true); try { @@ -309,55 +210,17 @@ FunctionParser::init_muparser() const // we may find after function names std::string transformed_expression = expressions[component]; - const char *function_names[] = {// functions predefined by muparser - "sin", - "cos", - "tan", - "asin", - "acos", - "atan", - "sinh", - "cosh", - "tanh", - "asinh", - "acosh", - "atanh", - "atan2", - "log2", - "log10", - "log", - "ln", - "exp", - "sqrt", - "sign", - "rint", - "abs", - "min", - "max", - "sum", - "avg", - // functions we define ourselves above - "if", - "int", - "ceil", - "cot", - "csc", - "floor", - "sec", - "pow", - "erfc", - "rand", - "rand_seed"}; - for (const auto &function_name_c_string : function_names) + for (const auto ¤t_function_name : + internal::FunctionParser::function_names) { - const std::string function_name = function_name_c_string; - const unsigned int function_name_length = function_name.size(); + const unsigned int function_name_length = + current_function_name.size(); std::string::size_type pos = 0; while (true) { // try to find any occurrences of the function name - pos = transformed_expression.find(function_name, pos); + pos = transformed_expression.find(current_function_name, pos); if (pos == std::string::npos) break; diff --git a/source/base/mu_parser_internal.cc b/source/base/mu_parser_internal.cc new file mode 100644 index 0000000000..9c4b964f23 --- /dev/null +++ b/source/base/mu_parser_internal.cc @@ -0,0 +1,203 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include + +#include + +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + + +#ifdef DEAL_II_WITH_MUPARSER + +namespace internal +{ + namespace FunctionParser + { + int + mu_round(double val) + { + return static_cast(val + ((val >= 0.0) ? 0.5 : -0.5)); + } + + double + mu_if(double condition, double thenvalue, double elsevalue) + { + if (mu_round(condition)) + return thenvalue; + else + return elsevalue; + } + + double + mu_or(double left, double right) + { + return (mu_round(left)) || (mu_round(right)); + } + + double + mu_and(double left, double right) + { + return (mu_round(left)) && (mu_round(right)); + } + + double + mu_int(double value) + { + return static_cast(mu_round(value)); + } + + double + mu_ceil(double value) + { + return std::ceil(value); + } + + double + mu_floor(double value) + { + return std::floor(value); + } + + double + mu_cot(double value) + { + return 1.0 / std::tan(value); + } + + double + mu_csc(double value) + { + return 1.0 / std::sin(value); + } + + double + mu_sec(double value) + { + return 1.0 / std::cos(value); + } + + double + mu_log(double value) + { + return std::log(value); + } + + double + mu_pow(double a, double b) + { + return std::pow(a, b); + } + + double + mu_erfc(double value) + { + return std::erfc(value); + } + + // returns a random value in the range [0,1] initializing the generator + // with the given seed + double + mu_rand_seed(double seed) + { + static Threads::Mutex rand_mutex; + std::lock_guard lock(rand_mutex); + + static boost::random::uniform_real_distribution<> uniform_distribution(0, + 1); + + // for each seed an unique random number generator is created, + // which is initialized with the seed itself + // we could use std::mt19937 but doing so results in compiler-dependent + // output. + static std::map rng_map; + + if (rng_map.find(seed) == rng_map.end()) + rng_map[seed] = boost::random::mt19937(static_cast(seed)); + + return uniform_distribution(rng_map[seed]); + } + + // returns a random value in the range [0,1] + double + mu_rand() + { + static Threads::Mutex rand_mutex; + std::lock_guard lock(rand_mutex); + static boost::random::uniform_real_distribution<> uniform_distribution(0, + 1); + static boost::random::mt19937 rng( + static_cast(std::time(nullptr))); + return uniform_distribution(rng); + } + + std::vector function_names = { + // functions predefined by muparser + "sin", + "cos", + "tan", + "asin", + "acos", + "atan", + "sinh", + "cosh", + "tanh", + "asinh", + "acosh", + "atanh", + "atan2", + "log2", + "log10", + "log", + "ln", + "exp", + "sqrt", + "sign", + "rint", + "abs", + "min", + "max", + "sum", + "avg", + // functions we define ourselves above + "if", + "int", + "ceil", + "cot", + "csc", + "floor", + "sec", + "pow", + "erfc", + "rand", + "rand_seed"}; + + } // namespace FunctionParser + +} // namespace internal +#endif + + + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/tensor_function_parser.cc b/source/base/tensor_function_parser.cc new file mode 100644 index 0000000000..df16de3589 --- /dev/null +++ b/source/base/tensor_function_parser.cc @@ -0,0 +1,400 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#ifdef DEAL_II_WITH_MUPARSER +# include +#endif + +DEAL_II_NAMESPACE_OPEN + + +template +const std::vector & +TensorFunctionParser::get_expressions() const +{ + return expressions; +} + + + +template +TensorFunctionParser::TensorFunctionParser( + const double initial_time) + : TensorFunction(initial_time) + , initialized(false) + , n_vars(0) + , n_components(Utilities::pow(dim, rank)) +{} + + +template +TensorFunctionParser::TensorFunctionParser( + const std::string &expression, + const std::string &constants, + const std::string &variable_names) + : TensorFunction() + , initialized(false) + , n_vars(0) + , n_components(Utilities::pow(dim, rank)) +{ + auto constants_map = Patterns::Tools::Convert::to_value( + constants, + std_cxx14::make_unique(Patterns::Anything(), + Patterns::Double(), + 0, + Patterns::Map::max_int_value, + ",", + "=")); + initialize(variable_names, + expression, + constants_map, + Utilities::split_string_list(variable_names, ",").size() == + dim + 1); +} + + +// We deliberately delay the definition of the default destructor +// so that we don't need to include the definition of mu::Parser +// in the header file. +template +TensorFunctionParser::~TensorFunctionParser() = default; + + +#ifdef DEAL_II_WITH_MUPARSER + +template +void +TensorFunctionParser::initialize( + const std::string & variables, + const std::vector & expressions, + const std::map &constants, + const bool time_dependent) +{ + this->tfp.clear(); // this will reset all thread-local objects + + this->constants = constants; + this->var_names = Utilities::split_string_list(variables, ','); + this->expressions = expressions; + AssertThrow(((time_dependent) ? dim + 1 : dim) == var_names.size(), + ExcMessage("Wrong number of variables")); + + // We check that the number of + // components of this function + // matches the number of components + // passed in as a vector of + // strings. + AssertThrow(this->n_components == expressions.size(), + ExcInvalidExpressionSize(this->n_components, expressions.size())); + + // Now we define how many variables + // we expect to read in. We + // distinguish between two cases: + // Time dependent problems, and not + // time dependent problems. In the + // first case the number of + // variables is given by the + // dimension plus one. In the other + // case, the number of variables is + // equal to the dimension. Once we + // parsed the variables string, if + // none of this is the case, then + // an exception is thrown. + if (time_dependent) + n_vars = dim + 1; + else + n_vars = dim; + + // create a parser object for the current thread we can then query + // in value() and vector_value(). this is not strictly necessary + // because a user may never call these functions on the current + // thread, but it gets us error messages about wrong formulas right + // away + init_muparser(); + + // finally set the initialization bit + initialized = true; +} + + + +template +void +TensorFunctionParser::init_muparser() const +{ + // check that we have not already initialized the parser on the + // current thread, i.e., that the current function is only called + // once per thread + Assert(tfp.get().size() == 0, ExcInternalError()); + + // initialize the objects for the current thread (tfp.get() and + // vars.get()) + tfp.get().reserve(this->n_components); + vars.get().resize(var_names.size()); + for (unsigned int component = 0; component < this->n_components; ++component) + { + tfp.get().emplace_back(new mu::Parser()); + + for (const auto &constant : constants) + { + tfp.get()[component]->DefineConst(constant.first, constant.second); + } + + for (unsigned int iv = 0; iv < var_names.size(); ++iv) + tfp.get()[component]->DefineVar(var_names[iv], &vars.get()[iv]); + + // define some compatibility functions: + tfp.get()[component]->DefineFun("if", + internal::FunctionParser::mu_if, + true); + tfp.get()[component]->DefineOprt("|", internal::FunctionParser::mu_or, 1); + tfp.get()[component]->DefineOprt("&", + internal::FunctionParser::mu_and, + 2); + tfp.get()[component]->DefineFun("int", + internal::FunctionParser::mu_int, + true); + tfp.get()[component]->DefineFun("ceil", + internal::FunctionParser::mu_ceil, + true); + tfp.get()[component]->DefineFun("cot", + internal::FunctionParser::mu_cot, + true); + tfp.get()[component]->DefineFun("csc", + internal::FunctionParser::mu_csc, + true); + tfp.get()[component]->DefineFun("floor", + internal::FunctionParser::mu_floor, + true); + tfp.get()[component]->DefineFun("sec", + internal::FunctionParser::mu_sec, + true); + tfp.get()[component]->DefineFun("log", + internal::FunctionParser::mu_log, + true); + tfp.get()[component]->DefineFun("pow", + internal::FunctionParser::mu_pow, + true); + tfp.get()[component]->DefineFun("erfc", + internal::FunctionParser::mu_erfc, + true); + tfp.get()[component]->DefineFun("rand_seed", + internal::FunctionParser::mu_rand_seed, + true); + tfp.get()[component]->DefineFun("rand", + internal::FunctionParser::mu_rand, + true); + + try + { + // muparser expects that functions have no + // space between the name of the function and the opening + // parenthesis. this is awkward because it is not backward + // compatible to the library we used to use before muparser + // (the tfparser library) but also makes no real sense. + // consequently, in the expressions we set, remove any space + // we may find after function names + std::string transformed_expression = expressions[component]; + + for (const auto ¤t_function_name : + internal::FunctionParser::function_names) + { + const unsigned int function_name_length = + current_function_name.size(); + + std::string::size_type pos = 0; + while (true) + { + // try to find any occurrences of the function name + pos = transformed_expression.find(current_function_name, pos); + if (pos == std::string::npos) + break; + + // replace whitespace until there no longer is any + while ((pos + function_name_length < + transformed_expression.size()) && + ((transformed_expression[pos + function_name_length] == + ' ') || + (transformed_expression[pos + function_name_length] == + '\t'))) + transformed_expression.erase( + transformed_expression.begin() + pos + + function_name_length); + + // move the current search position by the size of the + // actual function name + pos += function_name_length; + } + } + + // now use the transformed expression + tfp.get()[component]->SetExpr(transformed_expression); + } + catch (mu::ParserError &e) + { + std::cerr << "Message: <" << e.GetMsg() << ">\n"; + std::cerr << "Formula: <" << e.GetExpr() << ">\n"; + std::cerr << "Token: <" << e.GetToken() << ">\n"; + std::cerr << "Position: <" << e.GetPos() << ">\n"; + std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl; + AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg())); + } + } +} + + + +template +void +TensorFunctionParser::initialize( + const std::string & vars, + const std::string & expression, + const std::map &constants, + const bool time_dependent) +{ + initialize(vars, + Utilities::split_string_list(expression, ';'), + constants, + time_dependent); +} + + + +template +Tensor +TensorFunctionParser::value(const Point &p) const +{ + Assert(initialized == true, ExcNotInitialized()); + + // initialize the parser if that hasn't happened yet on the current thread + if (tfp.get().size() == 0) + init_muparser(); + + for (unsigned int i = 0; i < dim; ++i) + vars.get()[i] = p(i); + if (dim != n_vars) + vars.get()[dim] = this->get_time(); + + // initialize tensor with zeros + Tensor value; + + try + { + unsigned int component = 0; + for (Number *value_ptr = value.begin_raw(); value_ptr != value.end_raw(); + ++value_ptr) + { + *value_ptr = tfp.get()[component]->Eval(); + ++component; + } // for + } // try + catch (mu::ParserError &e) + { + std::cerr << "Message: <" << e.GetMsg() << ">\n"; + std::cerr << "Formula: <" << e.GetExpr() << ">\n"; + std::cerr << "Token: <" << e.GetToken() << ">\n"; + std::cerr << "Position: <" << e.GetPos() << ">\n"; + std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl; + AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg())); + + return Tensor(); + } // catch + + return value; +} + + + +template +void +TensorFunctionParser::value_list( + const std::vector> & p, + std::vector> &values) const +{ + Assert(p.size() == values.size(), + ExcDimensionMismatch(p.size(), values.size())); + + for (unsigned int i = 0; i < p.size(); ++i) + { + values[i] = value(p[i]); + } +} + +#else + + +template +void +TensorFunctionParser::initialize( + const std::string &, + const std::vector &, + const std::map &, + const bool) +{ + Assert(false, ExcNeedsFunctionparser()); +} + +template +void +TensorFunctionParser::initialize( + const std::string &, + const std::string &, + const std::map &, + const bool) +{ + Assert(false, ExcNeedsFunctionparser()); +} + + + +template +Tensor +TensorFunctionParser::value(const Point &) const +{ + Assert(false, ExcNeedsFunctionparser()); + return Tensor(); +} + + + +template +void +TensorFunctionParser::value_list( + const std::vector> & p, + std::vector> &values) const +{ + Assert(false, ExcNeedsFunctionparser()); +} + + +#endif + +// explicit instantiations +#include "tensor_function_parser.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/tensor_function_parser.inst.in b/source/base/tensor_function_parser.inst.in new file mode 100644 index 0000000000..3e73de5f02 --- /dev/null +++ b/source/base/tensor_function_parser.inst.in @@ -0,0 +1,25 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +for (S : REAL_SCALARS; rank : RANKS; dim : SPACE_DIMENSIONS) + { + template class TensorFunctionParser; + } + +for (S : COMPLEX_SCALARS; rank : RANKS; dim : SPACE_DIMENSIONS) + { + template class TensorFunctionParser; + } diff --git a/tests/base/tensor_function_parser_01.cc b/tests/base/tensor_function_parser_01.cc new file mode 100644 index 0000000000..c311c865fd --- /dev/null +++ b/tests/base/tensor_function_parser_01.cc @@ -0,0 +1,69 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// This program tests the functionality of the function parser +// wrapper. + +#include +#include +#include + +#include "../tests.h" + + +// Test initializtaion and evaluation of timedependent function in 2d. +int +main() +{ + initlog(); + + // Define some constants that will be used by the function parser + std::map constants; + constants["pi"] = numbers::PI; + + // Define the variables that will be used inside the expressions + std::string variables = "x,y,t"; + + // Define the expressions of the vector_valued function. + std::vector expressions; + expressions.push_back("cos(2*pi*(x*y+t))"); + expressions.push_back("sin(2*pi*(x*y+t))"); + expressions.push_back("-sin(2*pi*(x*y+t))"); + expressions.push_back("cos(2*pi*(x*y+t))"); + + // Test time dependent function + try + { + { + TensorFunctionParser<2, 2> tensor_function; + tensor_function.initialize(variables, + expressions, + constants, + /* time dependent */ true); + + deallog << "Initialize Succeeded with dim = 2, rank = 2, " + << expressions.size() << " expressions, " << variables + << " as variables." << std::endl; + } + } + catch (...) + { + deallog << "Initialization or Evaluation Failed with dim = 2, rank = 2, " + << expressions.size() << " expressions, " << variables + << " as variables." << std::endl; + } +} diff --git a/tests/base/tensor_function_parser_01.with_muparser=true.output b/tests/base/tensor_function_parser_01.with_muparser=true.output new file mode 100644 index 0000000000..57b39e740b --- /dev/null +++ b/tests/base/tensor_function_parser_01.with_muparser=true.output @@ -0,0 +1,2 @@ + +DEAL::Initialize Succeeded with dim = 2, rank = 2, 4 expressions, x,y,t as variables. diff --git a/tests/base/tensor_function_parser_02.cc b/tests/base/tensor_function_parser_02.cc new file mode 100644 index 0000000000..4081ba8d7f --- /dev/null +++ b/tests/base/tensor_function_parser_02.cc @@ -0,0 +1,69 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// This program tests the functionality of the function parser +// wrapper. + +#include +#include +#include + +#include "../tests.h" + + +// Test initializtaion and evaluation of timedependent function in 2d. +int +main() +{ + initlog(); + + // Define some constants that will be used by the function parser + std::map constants; + constants["pi"] = numbers::PI; + + // Define the variables that will be used inside the expressions + std::string variables = "x,y"; + + // Define the expressions of the vector_valued function. + std::vector expressions; + expressions.push_back("cos(2*pi*(x*y))"); + expressions.push_back("sin(2*pi*(x*y))"); + expressions.push_back("-sin(2*pi*(x*y))"); + expressions.push_back("cos(2*pi*(x*y))"); + + // Test time dependent function + try + { + { + TensorFunctionParser<2, 2> tensor_function; + tensor_function.initialize(variables, + expressions, + constants, + /* time dependent */ false); + + deallog << "Initialize Succeeded with dim = 2, rank = 2, " + << expressions.size() << " expressions, " << variables + << " as variables." << std::endl; + } + } + catch (...) + { + deallog << "Initialization or Evaluation Failed with dim = 2, rank = 2, " + << expressions.size() << " expressions, " << variables + << " as variables." << std::endl; + } +} diff --git a/tests/base/tensor_function_parser_02.with_muparser=true.output b/tests/base/tensor_function_parser_02.with_muparser=true.output new file mode 100644 index 0000000000..5fb16609c6 --- /dev/null +++ b/tests/base/tensor_function_parser_02.with_muparser=true.output @@ -0,0 +1,2 @@ + +DEAL::Initialize Succeeded with dim = 2, rank = 2, 4 expressions, x,y as variables. diff --git a/tests/base/tensor_function_parser_03.cc b/tests/base/tensor_function_parser_03.cc new file mode 100644 index 0000000000..7bcc249c9c --- /dev/null +++ b/tests/base/tensor_function_parser_03.cc @@ -0,0 +1,82 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// This program tests the functionality of the function parser +// wrapper. + +#include +#include +#include + +#include "../tests.h" + + +// Test initializtaion and evaluation of timedependent function in 2d. +int +main() +{ + initlog(); + + // Define some constants that will be used by the function parser + std::map constants; + constants["pi"] = numbers::PI; + + // Define the variables that will be used inside the expressions + std::string variables = "x,y,t"; + + // Define the expressions of the vector_valued function. + std::vector expressions; + expressions.push_back("cos(2*pi*(x*y+t))"); + expressions.push_back("sin(2*pi*(x*y+t))"); + expressions.push_back("-sin(2*pi*(x*y+t))"); + expressions.push_back("cos(2*pi*(x*y+t))"); + + Point<2> p(0.5, 0.5); // evaluation point + + // Test time dependent function + try + { + { + TensorFunctionParser<2, 2> tensor_function; + tensor_function.initialize(variables, + expressions, + constants, + /* time dependent */ true); + + deallog << "Initialize Succeeded with dim = 2, rank = 2, " + << expressions.size() << " expressions, " << variables + << " as variables." << std::endl; + + for (unsigned int t = 0; t < 5; ++t) + { + tensor_function.set_time(t); + deallog << "Value at p = (0.5,0.5): " << std::endl + << tensor_function.value(p)[0][0] << " " + << tensor_function.value(p)[0][1] << std::endl + << tensor_function.value(p)[1][0] << " " + << tensor_function.value(p)[1][1] << " at time t = " << t + << std::endl; + } + } + } + catch (...) + { + deallog << "Initialization or Evaluation Failed with dim = 2, rank = 2, " + << expressions.size() << " expressions, " << variables + << " as variables." << std::endl; + } +} diff --git a/tests/base/tensor_function_parser_03.with_muparser=true.output b/tests/base/tensor_function_parser_03.with_muparser=true.output new file mode 100644 index 0000000000..7a64cc8a28 --- /dev/null +++ b/tests/base/tensor_function_parser_03.with_muparser=true.output @@ -0,0 +1,17 @@ + +DEAL::Initialize Succeeded with dim = 2, rank = 2, 4 expressions, x,y,t as variables. +DEAL::Value at p = (0.5,0.5): +DEAL::6.12323e-17 1.00000 +DEAL::-1.00000 6.12323e-17 at time t = 0 +DEAL::Value at p = (0.5,0.5): +DEAL::3.06162e-16 1.00000 +DEAL::-1.00000 3.06162e-16 at time t = 1 +DEAL::Value at p = (0.5,0.5): +DEAL::5.51091e-16 1.00000 +DEAL::-1.00000 5.51091e-16 at time t = 2 +DEAL::Value at p = (0.5,0.5): +DEAL::-9.80336e-16 1.00000 +DEAL::-1.00000 -9.80336e-16 at time t = 3 +DEAL::Value at p = (0.5,0.5): +DEAL::-7.35407e-16 1.00000 +DEAL::-1.00000 -7.35407e-16 at time t = 4 diff --git a/tests/base/tensor_function_parser_04.cc b/tests/base/tensor_function_parser_04.cc new file mode 100644 index 0000000000..2f0fd76b6b --- /dev/null +++ b/tests/base/tensor_function_parser_04.cc @@ -0,0 +1,78 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// This program tests the functionality of the function parser +// wrapper. + +#include +#include +#include + +#include "../tests.h" + + +// Test initializtaion and evaluation of timedependent function in 2d. +int +main() +{ + initlog(); + + // Define some constants that will be used by the function parser + std::map constants; + constants["pi"] = numbers::PI; + + // Define the variables that will be used inside the expressions + std::string variables = "x,y,t"; + + // Define the expressions of the vector_valued function. + std::vector expressions; + expressions.push_back("cos(2*pi*(x*y+t+0.5))"); + expressions.push_back("sin(2*pi*(x*y-t+0.5))"); + + Point<2> p(0.5, 0.5); // evaluation point + + // Test time dependent function + try + { + { + TensorFunctionParser<1, 2> tensor_function; + tensor_function.initialize(variables, + expressions, + constants, + /* time dependent */ true); + + deallog << "Initialize Succeeded with dim = 2, rank = 1, " + << expressions.size() << " expressions, " << variables + << " as variables." << std::endl; + + for (unsigned int t = 0; t < 5; ++t) + { + tensor_function.set_time(t); + deallog << "Value at p = (0.5,0.5): " << std::endl + << tensor_function.value(p)[0] << " " + << tensor_function.value(p)[1] << " at time t = " << t + << std::endl; + } + } + } + catch (...) + { + deallog << "Initialization or Evaluation Failed with dim = 2, rank = 2, " + << expressions.size() << " expressions, " << variables + << " as variables." << std::endl; + } +} diff --git a/tests/base/tensor_function_parser_04.with_muparser=true.output b/tests/base/tensor_function_parser_04.with_muparser=true.output new file mode 100644 index 0000000000..a81d0acbbe --- /dev/null +++ b/tests/base/tensor_function_parser_04.with_muparser=true.output @@ -0,0 +1,12 @@ + +DEAL::Initialize Succeeded with dim = 2, rank = 1, 2 expressions, x,y,t as variables. +DEAL::Value at p = (0.5,0.5): +DEAL::-1.83697e-16 -1.00000 at time t = 0 +DEAL::Value at p = (0.5,0.5): +DEAL::-4.28626e-16 -1.00000 at time t = 1 +DEAL::Value at p = (0.5,0.5): +DEAL::-2.44991e-15 -1.00000 at time t = 2 +DEAL::Value at p = (0.5,0.5): +DEAL::-2.69484e-15 -1.00000 at time t = 3 +DEAL::Value at p = (0.5,0.5): +DEAL::-2.93977e-15 -1.00000 at time t = 4 -- 2.39.5