Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
function_parser.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/function_parser.h>
18 #include <deal.II/base/patterns.h>
19 #include <deal.II/base/thread_management.h>
20 #include <deal.II/base/utilities.h>
21 
22 #include <deal.II/lac/vector.h>
23 
24 #include <boost/random.hpp>
25 
26 #include <cmath>
27 #include <map>
28 
29 #ifdef DEAL_II_WITH_MUPARSER
30 # include <muParser.h>
31 #else
32 
33 
34 
35 namespace fparser
36 {
37  class FunctionParser
38  {};
39 } // namespace fparser
40 #endif
41 
42 DEAL_II_NAMESPACE_OPEN
43 
44 
45 template <int dim>
46 const std::vector<std::string> &
48 {
49  return expressions;
50 }
51 
52 
53 
54 template <int dim>
56  const double initial_time,
57  const double h)
58  : AutoDerivativeFunction<dim>(h, n_components, initial_time)
59  , initialized(false)
60  , n_vars(0)
61 {}
62 
63 
64 template <int dim>
65 FunctionParser<dim>::FunctionParser(const std::string &expression,
66  const std::string &constants,
67  const std::string &variable_names,
68  const double h)
70  h,
71  Utilities::split_string_list(expression, ';').size())
72  , initialized(false)
73  , n_vars(0)
74 {
76  constants,
77  std_cxx14::make_unique<Patterns::Map>(Patterns::Anything(),
79  0,
81  ",",
82  "="));
83  initialize(variable_names,
84  expression,
85  constants_map,
86  Utilities::split_string_list(variable_names, ",").size() ==
87  dim + 1);
88 }
89 
90 
91 // We deliberately delay the definition of the default destructor
92 // so that we don't need to include the definition of mu::Parser
93 // in the header file.
94 template <int dim>
96 
97 
98 #ifdef DEAL_II_WITH_MUPARSER
99 
100 template <int dim>
101 void
102 FunctionParser<dim>::initialize(const std::string & variables,
103  const std::vector<std::string> &expressions,
104  const std::map<std::string, double> &constants,
105  const bool time_dependent)
106 {
107  this->fp.clear(); // this will reset all thread-local objects
108 
109  this->constants = constants;
110  this->var_names = Utilities::split_string_list(variables, ',');
111  this->expressions = expressions;
112  AssertThrow(((time_dependent) ? dim + 1 : dim) == var_names.size(),
113  ExcMessage("Wrong number of variables"));
114 
115  // We check that the number of
116  // components of this function
117  // matches the number of components
118  // passed in as a vector of
119  // strings.
120  AssertThrow(this->n_components == expressions.size(),
121  ExcInvalidExpressionSize(this->n_components, expressions.size()));
122 
123  // Now we define how many variables
124  // we expect to read in. We
125  // distinguish between two cases:
126  // Time dependent problems, and not
127  // time dependent problems. In the
128  // first case the number of
129  // variables is given by the
130  // dimension plus one. In the other
131  // case, the number of variables is
132  // equal to the dimension. Once we
133  // parsed the variables string, if
134  // none of this is the case, then
135  // an exception is thrown.
136  if (time_dependent)
137  n_vars = dim + 1;
138  else
139  n_vars = dim;
140 
141  // create a parser object for the current thread we can then query
142  // in value() and vector_value(). this is not strictly necessary
143  // because a user may never call these functions on the current
144  // thread, but it gets us error messages about wrong formulas right
145  // away
146  init_muparser();
147 
148  // finally set the initialization bit
149  initialized = true;
150 }
151 
152 
153 
154 namespace internal
155 {
156  // convert double into int
157  int
158  mu_round(double val)
159  {
160  return static_cast<int>(val + ((val >= 0.0) ? 0.5 : -0.5));
161  }
162 
163  double
164  mu_if(double condition, double thenvalue, double elsevalue)
165  {
166  if (mu_round(condition))
167  return thenvalue;
168  else
169  return elsevalue;
170  }
171 
172  double
173  mu_or(double left, double right)
174  {
175  return (mu_round(left)) || (mu_round(right));
176  }
177 
178  double
179  mu_and(double left, double right)
180  {
181  return (mu_round(left)) && (mu_round(right));
182  }
183 
184  double
185  mu_int(double value)
186  {
187  return static_cast<double>(mu_round(value));
188  }
189 
190  double
191  mu_ceil(double value)
192  {
193  return std::ceil(value);
194  }
195 
196  double
197  mu_floor(double value)
198  {
199  return std::floor(value);
200  }
201 
202  double
203  mu_cot(double value)
204  {
205  return 1.0 / std::tan(value);
206  }
207 
208  double
209  mu_csc(double value)
210  {
211  return 1.0 / std::sin(value);
212  }
213 
214  double
215  mu_sec(double value)
216  {
217  return 1.0 / std::cos(value);
218  }
219 
220  double
221  mu_log(double value)
222  {
223  return std::log(value);
224  }
225 
226  double
227  mu_pow(double a, double b)
228  {
229  return std::pow(a, b);
230  }
231 
232  double
233  mu_erfc(double value)
234  {
235  return std::erfc(value);
236  }
237 
238  // returns a random value in the range [0,1] initializing the generator
239  // with the given seed
240  double
241  mu_rand_seed(double seed)
242  {
243  static Threads::Mutex rand_mutex;
244  std::lock_guard<std::mutex> lock(rand_mutex);
245 
246  static boost::random::uniform_real_distribution<> uniform_distribution(0,
247  1);
248 
249  // for each seed an unique random number generator is created,
250  // which is initialized with the seed itself
251  static std::map<double, boost::random::mt19937> rng_map;
252 
253  if (rng_map.find(seed) == rng_map.end())
254  rng_map[seed] = boost::random::mt19937(static_cast<unsigned int>(seed));
255 
256  return uniform_distribution(rng_map[seed]);
257  }
258 
259  // returns a random value in the range [0,1]
260  double
261  mu_rand()
262  {
263  static Threads::Mutex rand_mutex;
264  std::lock_guard<std::mutex> lock(rand_mutex);
265  static boost::random::uniform_real_distribution<> uniform_distribution(0,
266  1);
267  static boost::random::mt19937 rng(
268  static_cast<unsigned long>(std::time(nullptr)));
269  return uniform_distribution(rng);
270  }
271 
272 } // namespace internal
273 
274 
275 
276 template <int dim>
277 void
279 {
280  // check that we have not already initialized the parser on the
281  // current thread, i.e., that the current function is only called
282  // once per thread
283  Assert(fp.get().size() == 0, ExcInternalError());
284 
285  // initialize the objects for the current thread (fp.get() and
286  // vars.get())
287  fp.get().reserve(this->n_components);
288  vars.get().resize(var_names.size());
289  for (unsigned int component = 0; component < this->n_components; ++component)
290  {
291  fp.get().emplace_back(new mu::Parser());
292 
293  for (std::map<std::string, double>::const_iterator constant =
294  constants.begin();
295  constant != constants.end();
296  ++constant)
297  {
298  fp.get()[component]->DefineConst(constant->first, constant->second);
299  }
300 
301  for (unsigned int iv = 0; iv < var_names.size(); ++iv)
302  fp.get()[component]->DefineVar(var_names[iv], &vars.get()[iv]);
303 
304  // define some compatibility functions:
305  fp.get()[component]->DefineFun("if", internal::mu_if, true);
306  fp.get()[component]->DefineOprt("|", internal::mu_or, 1);
307  fp.get()[component]->DefineOprt("&", internal::mu_and, 2);
308  fp.get()[component]->DefineFun("int", internal::mu_int, true);
309  fp.get()[component]->DefineFun("ceil", internal::mu_ceil, true);
310  fp.get()[component]->DefineFun("cot", internal::mu_cot, true);
311  fp.get()[component]->DefineFun("csc", internal::mu_csc, true);
312  fp.get()[component]->DefineFun("floor", internal::mu_floor, true);
313  fp.get()[component]->DefineFun("sec", internal::mu_sec, true);
314  fp.get()[component]->DefineFun("log", internal::mu_log, true);
315  fp.get()[component]->DefineFun("pow", internal::mu_pow, true);
316  fp.get()[component]->DefineFun("erfc", internal::mu_erfc, true);
317  fp.get()[component]->DefineFun("rand_seed", internal::mu_rand_seed, true);
318  fp.get()[component]->DefineFun("rand", internal::mu_rand, true);
319 
320  try
321  {
322  // muparser expects that functions have no
323  // space between the name of the function and the opening
324  // parenthesis. this is awkward because it is not backward
325  // compatible to the library we used to use before muparser
326  // (the fparser library) but also makes no real sense.
327  // consequently, in the expressions we set, remove any space
328  // we may find after function names
329  std::string transformed_expression = expressions[component];
330 
331  const char *function_names[] = {// functions predefined by muparser
332  "sin",
333  "cos",
334  "tan",
335  "asin",
336  "acos",
337  "atan",
338  "sinh",
339  "cosh",
340  "tanh",
341  "asinh",
342  "acosh",
343  "atanh",
344  "atan2",
345  "log2",
346  "log10",
347  "log",
348  "ln",
349  "exp",
350  "sqrt",
351  "sign",
352  "rint",
353  "abs",
354  "min",
355  "max",
356  "sum",
357  "avg",
358  // functions we define ourselves above
359  "if",
360  "int",
361  "ceil",
362  "cot",
363  "csc",
364  "floor",
365  "sec",
366  "pow",
367  "erfc",
368  "rand",
369  "rand_seed"};
370  for (const auto &function_name_c_string : function_names)
371  {
372  const std::string function_name = function_name_c_string;
373  const unsigned int function_name_length = function_name.size();
374 
375  std::string::size_type pos = 0;
376  while (true)
377  {
378  // try to find any occurrences of the function name
379  pos = transformed_expression.find(function_name, pos);
380  if (pos == std::string::npos)
381  break;
382 
383  // replace whitespace until there no longer is any
384  while ((pos + function_name_length <
385  transformed_expression.size()) &&
386  ((transformed_expression[pos + function_name_length] ==
387  ' ') ||
388  (transformed_expression[pos + function_name_length] ==
389  '\t')))
390  transformed_expression.erase(
391  transformed_expression.begin() + pos +
392  function_name_length);
393 
394  // move the current search position by the size of the
395  // actual function name
396  pos += function_name_length;
397  }
398  }
399 
400  // now use the transformed expression
401  fp.get()[component]->SetExpr(transformed_expression);
402  }
403  catch (mu::ParserError &e)
404  {
405  std::cerr << "Message: <" << e.GetMsg() << ">\n";
406  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
407  std::cerr << "Token: <" << e.GetToken() << ">\n";
408  std::cerr << "Position: <" << e.GetPos() << ">\n";
409  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
410  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
411  }
412  }
413 }
414 
415 
416 
417 template <int dim>
418 void
419 FunctionParser<dim>::initialize(const std::string & vars,
420  const std::string & expression,
421  const std::map<std::string, double> &constants,
422  const bool time_dependent)
423 {
424  initialize(vars,
425  Utilities::split_string_list(expression, ';'),
426  constants,
427  time_dependent);
428 }
429 
430 
431 
432 template <int dim>
433 double
435  const unsigned int component) const
436 {
437  Assert(initialized == true, ExcNotInitialized());
438  Assert(component < this->n_components,
439  ExcIndexRange(component, 0, this->n_components));
440 
441  // initialize the parser if that hasn't happened yet on the current thread
442  if (fp.get().size() == 0)
443  init_muparser();
444 
445  for (unsigned int i = 0; i < dim; ++i)
446  vars.get()[i] = p(i);
447  if (dim != n_vars)
448  vars.get()[dim] = this->get_time();
449 
450  try
451  {
452  return fp.get()[component]->Eval();
453  }
454  catch (mu::ParserError &e)
455  {
456  std::cerr << "Message: <" << e.GetMsg() << ">\n";
457  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
458  std::cerr << "Token: <" << e.GetToken() << ">\n";
459  std::cerr << "Position: <" << e.GetPos() << ">\n";
460  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
461  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
462  return 0.0;
463  }
464 }
465 
466 
467 
468 template <int dim>
469 void
471  Vector<double> & values) const
472 {
473  Assert(initialized == true, ExcNotInitialized());
474  Assert(values.size() == this->n_components,
475  ExcDimensionMismatch(values.size(), this->n_components));
476 
477 
478  // initialize the parser if that hasn't happened yet on the current thread
479  if (fp.get().size() == 0)
480  init_muparser();
481 
482  for (unsigned int i = 0; i < dim; ++i)
483  vars.get()[i] = p(i);
484  if (dim != n_vars)
485  vars.get()[dim] = this->get_time();
486 
487  for (unsigned int component = 0; component < this->n_components; ++component)
488  values(component) = fp.get()[component]->Eval();
489 }
490 
491 #else
492 
493 
494 template <int dim>
495 void
496 FunctionParser<dim>::initialize(const std::string &,
497  const std::vector<std::string> &,
498  const std::map<std::string, double> &,
499  const bool)
500 {
501  Assert(false, ExcNeedsFunctionparser());
502 }
503 
504 template <int dim>
505 void
506 FunctionParser<dim>::initialize(const std::string &,
507  const std::string &,
508  const std::map<std::string, double> &,
509  const bool)
510 {
511  Assert(false, ExcNeedsFunctionparser());
512 }
513 
514 
515 
516 template <int dim>
517 double
518 FunctionParser<dim>::value(const Point<dim> &, unsigned int) const
519 {
520  Assert(false, ExcNeedsFunctionparser());
521  return 0.;
522 }
523 
524 
525 template <int dim>
526 void
528 {
529  Assert(false, ExcNeedsFunctionparser());
530 }
531 
532 
533 #endif
534 
535 // Explicit Instantiations.
536 
537 template class FunctionParser<1>;
538 template class FunctionParser<2>;
539 template class FunctionParser<3>;
540 
541 DEAL_II_NAMESPACE_CLOSE
std::vector< std::string > split_string_list(const std::string &s, const std::string &delimiter=",")
Definition: utilities.cc:566
static ::ExceptionBase & ExcNeedsFunctionparser()
FunctionParser(const unsigned int n_components=1, const double initial_time=0.0, const double h=1e-8)
static T to_value(const std::string &s, const std::unique_ptr< Patterns::PatternBase > &p=Convert< T >::to_pattern())=delete
static const unsigned int max_int_value
Definition: patterns.h:591
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void init_muparser() const
Definition: cuda.h:31
const std::vector< std::string > & get_expressions() const
void initialize(const std::string &vars, const std::vector< std::string > &expressions, const ConstMap &constants, const bool time_dependent=false)
size_type size() const
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
~FunctionParser() override
std::map< std::string, double > constants
static ::ExceptionBase & ExcInternalError()