Reference documentation for deal.II version GIT 1e9e64a91d 2022-09-28 19:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mu_parser_internal.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2022 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 
18 #include <deal.II/base/utilities.h>
19 
20 #include <cmath>
21 #include <ctime>
22 #include <limits>
23 #include <map>
24 #include <mutex>
25 #include <random>
26 #include <vector>
27 
28 #ifdef DEAL_II_WITH_MUPARSER
29 # include <muParser.h>
30 #endif
31 
33 
34 namespace internal
35 {
36  namespace FunctionParser
37  {
38  int
39  mu_round(double val)
40  {
41  return static_cast<int>(val + ((val >= 0.0) ? 0.5 : -0.5));
42  }
43 
44  double
45  mu_if(double condition, double thenvalue, double elsevalue)
46  {
47  if (mu_round(condition) != 0)
48  return thenvalue;
49  else
50  return elsevalue;
51  }
52 
53  double
54  mu_or(double left, double right)
55  {
56  return static_cast<double>((mu_round(left) != 0) ||
57  (mu_round(right) != 0));
58  }
59 
60  double
61  mu_and(double left, double right)
62  {
63  return static_cast<double>((mu_round(left) != 0) &&
64  (mu_round(right) != 0));
65  }
66 
67  double
68  mu_int(double value)
69  {
70  return static_cast<double>(mu_round(value));
71  }
72 
73  double
74  mu_ceil(double value)
75  {
76  return std::ceil(value);
77  }
78 
79  double
80  mu_floor(double value)
81  {
82  return std::floor(value);
83  }
84 
85  double
86  mu_cot(double value)
87  {
88  return 1.0 / std::tan(value);
89  }
90 
91  double
92  mu_csc(double value)
93  {
94  return 1.0 / std::sin(value);
95  }
96 
97  double
98  mu_sec(double value)
99  {
100  return 1.0 / std::cos(value);
101  }
102 
103  double
104  mu_log(double value)
105  {
106  return std::log(value);
107  }
108 
109  double
110  mu_pow(double a, double b)
111  {
112  return std::pow(a, b);
113  }
114 
115  double
116  mu_erf(double value)
117  {
118  return std::erf(value);
119  }
120 
121  double
122  mu_erfc(double value)
123  {
124  return std::erfc(value);
125  }
126 
127  // returns a random value in the range [0,1] initializing the generator
128  // with the given seed
129  double
130  mu_rand_seed(double seed)
131  {
132  static std::mutex rand_mutex;
133  std::lock_guard<std::mutex> lock(rand_mutex);
134 
135  std::uniform_real_distribution<> uniform_distribution(0., 1.);
136 
137  // for each seed a unique random number generator is created,
138  // which is initialized with the seed itself
139  static std::map<double, std::mt19937> rng_map;
140 
141  if (rng_map.find(seed) == rng_map.end())
142  rng_map[seed] = std::mt19937(static_cast<unsigned int>(seed));
143 
144  return uniform_distribution(rng_map[seed]);
145  }
146 
147  // returns a random value in the range [0,1]
148  double
150  {
151  static std::mutex rand_mutex;
152  std::lock_guard<std::mutex> lock(rand_mutex);
153  std::uniform_real_distribution<> uniform_distribution(0., 1.);
154  const unsigned int seed = static_cast<unsigned long>(std::time(nullptr));
155  static std::mt19937 rng(seed);
156  return uniform_distribution(rng);
157  }
158 
159  std::vector<std::string>
161  {
162  return {// functions predefined by muparser
163  "sin",
164  "cos",
165  "tan",
166  "asin",
167  "acos",
168  "atan",
169  "sinh",
170  "cosh",
171  "tanh",
172  "asinh",
173  "acosh",
174  "atanh",
175  "atan2",
176  "log2",
177  "log10",
178  "log",
179  "ln",
180  "exp",
181  "sqrt",
182  "sign",
183  "rint",
184  "abs",
185  "min",
186  "max",
187  "sum",
188  "avg",
189  // functions we define ourselves above
190  "if",
191  "int",
192  "ceil",
193  "cot",
194  "csc",
195  "floor",
196  "sec",
197  "pow",
198  "erf",
199  "erfc",
200  "rand",
201  "rand_seed"};
202  }
203 
204 #ifdef DEAL_II_WITH_MUPARSER
208  class Parser : public muParserBase
209  {
210  public:
211  operator mu::Parser &()
212  {
213  return parser;
214  }
215 
216  operator const mu::Parser &() const
217  {
218  return parser;
219  }
220 
221  protected:
222  mu::Parser parser;
223  };
224 #endif
225 
226  template <int dim, typename Number>
228  : initialized(false)
229  , n_vars(0)
230  {}
231 
232 
233 
234  template <int dim, typename Number>
235  void
237  const std::string & variables,
238  const std::vector<std::string> & expressions,
239  const std::map<std::string, double> &constants,
240  const bool time_dependent)
241  {
242  this->parser_data.clear(); // this will reset all thread-local objects
243 
244  this->constants = constants;
245  this->var_names = Utilities::split_string_list(variables, ',');
246  this->expressions = expressions;
247  AssertThrow(((time_dependent) ? dim + 1 : dim) == this->var_names.size(),
248  ExcMessage("Wrong number of variables"));
249 
250  // Now we define how many variables we expect to read in. We distinguish
251  // between two cases: Time dependent problems, and not time dependent
252  // problems. In the first case the number of variables is given by the
253  // dimension plus one. In the other case, the number of variables is equal
254  // to the dimension. Once we parsed the variables string, if none of this
255  // is the case, then an exception is thrown.
256  if (time_dependent)
257  this->n_vars = dim + 1;
258  else
259  this->n_vars = dim;
260 
261  // create a parser object for the current thread we can then query in
262  // value() and vector_value(). this is not strictly necessary because a
263  // user may never call these functions on the current thread, but it gets
264  // us error messages about wrong formulas right away
265  this->init_muparser();
266  this->initialized = true;
267  }
268 
269 
270 
271  template <int dim, typename Number>
272  void
274  {
275 #ifdef DEAL_II_WITH_MUPARSER
276  // check that we have not already initialized the parser on the
277  // current thread, i.e., that the current function is only called
278  // once per thread
279  ParserData &data = this->parser_data.get();
280  Assert(data.parsers.size() == 0 && data.vars.size() == 0,
281  ExcInternalError());
282  const unsigned int n_components = expressions.size();
283 
284  // initialize the objects for the current thread
285  data.parsers.reserve(n_components);
286  data.vars.resize(this->var_names.size());
287  for (unsigned int component = 0; component < n_components; ++component)
288  {
289  data.parsers.emplace_back(std::make_unique<Parser>());
290  mu::Parser &parser = dynamic_cast<Parser &>(*data.parsers.back());
291 
292  for (const auto &constant : this->constants)
293  parser.DefineConst(constant.first, constant.second);
294 
295  for (unsigned int iv = 0; iv < this->var_names.size(); ++iv)
296  parser.DefineVar(this->var_names[iv], &data.vars[iv]);
297 
298  // define some compatibility functions:
299  parser.DefineFun("if", mu_if, true);
300  parser.DefineOprt("|", mu_or, 1);
301  parser.DefineOprt("&", mu_and, 2);
302  parser.DefineFun("int", mu_int, true);
303  parser.DefineFun("ceil", mu_ceil, true);
304  parser.DefineFun("cot", mu_cot, true);
305  parser.DefineFun("csc", mu_csc, true);
306  parser.DefineFun("floor", mu_floor, true);
307  parser.DefineFun("sec", mu_sec, true);
308  parser.DefineFun("log", mu_log, true);
309  parser.DefineFun("pow", mu_pow, true);
310  parser.DefineFun("erfc", mu_erfc, true);
311  parser.DefineFun("rand_seed", mu_rand_seed, true);
312  parser.DefineFun("rand", mu_rand, true);
313 
314  try
315  {
316  // muparser expects that functions have no
317  // space between the name of the function and the opening
318  // parenthesis. this is awkward because it is not backward
319  // compatible to the library we used to use before muparser
320  // (the fparser library) but also makes no real sense.
321  // consequently, in the expressions we set, remove any space
322  // we may find after function names
323  std::string transformed_expression = this->expressions[component];
324 
325  for (const auto &current_function_name : get_function_names())
326  {
327  const unsigned int function_name_length =
328  current_function_name.size();
329 
330  std::string::size_type pos = 0;
331  while (true)
332  {
333  // try to find any occurrences of the function name
334  pos =
335  transformed_expression.find(current_function_name, pos);
336  if (pos == std::string::npos)
337  break;
338 
339  // replace whitespace until there no longer is any
340  while (
341  (pos + function_name_length <
342  transformed_expression.size()) &&
343  ((transformed_expression[pos + function_name_length] ==
344  ' ') ||
345  (transformed_expression[pos + function_name_length] ==
346  '\t')))
347  transformed_expression.erase(
348  transformed_expression.begin() + pos +
349  function_name_length);
350 
351  // move the current search position by the size of the
352  // actual function name
353  pos += function_name_length;
354  }
355  }
356 
357  // now use the transformed expression
358  parser.SetExpr(transformed_expression);
359  }
360  catch (mu::ParserError &e)
361  {
362  std::cerr << "Message: <" << e.GetMsg() << ">\n";
363  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
364  std::cerr << "Token: <" << e.GetToken() << ">\n";
365  std::cerr << "Position: <" << e.GetPos() << ">\n";
366  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
367  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
368  }
369  }
370 #else
372 #endif
373  }
374 
375  template <int dim, typename Number>
376  Number
378  const double time,
379  unsigned int component) const
380  {
381 #ifdef DEAL_II_WITH_MUPARSER
382  Assert(this->initialized == true, ExcNotInitialized());
383 
384  // initialize the parser if that hasn't happened yet on the current
385  // thread
386  internal::FunctionParser::ParserData &data = this->parser_data.get();
387  if (data.vars.size() == 0)
388  init_muparser();
389 
390  for (unsigned int i = 0; i < dim; ++i)
391  data.vars[i] = p(i);
392  if (dim != this->n_vars)
393  data.vars[dim] = time;
394 
395  try
396  {
397  Assert(dynamic_cast<Parser *>(data.parsers[component].get()),
398  ExcInternalError());
399  // NOLINTNEXTLINE don't warn about using static_cast once we check
400  mu::Parser &parser = static_cast<Parser &>(*data.parsers[component]);
401  return parser.Eval();
402  } // try
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  } // catch
412 
413 #else
414  (void)p;
415  (void)time;
416  (void)component;
418 #endif
419  return std::numeric_limits<double>::signaling_NaN();
420  }
421 
422  template <int dim, typename Number>
423  void
425  const Point<dim> & p,
426  const double time,
427  ArrayView<Number> &values) const
428  {
429 #ifdef DEAL_II_WITH_MUPARSER
430  Assert(this->initialized == true, ExcNotInitialized());
431 
432  // initialize the parser if that hasn't happened yet on the current
433  // thread
434  internal::FunctionParser::ParserData &data = this->parser_data.get();
435  if (data.vars.size() == 0)
436  init_muparser();
437 
438  for (unsigned int i = 0; i < dim; ++i)
439  data.vars[i] = p(i);
440  if (dim != this->n_vars)
441  data.vars[dim] = time;
442 
443  AssertDimension(values.size(), data.parsers.size());
444  try
445  {
446  for (unsigned int component = 0; component < data.parsers.size();
447  ++component)
448  {
449  Assert(dynamic_cast<Parser *>(data.parsers[component].get()),
450  ExcInternalError());
451  mu::Parser &parser =
452  // We just checked that the pointer is valid so suppress the
453  // clang-tidy check
454  static_cast<Parser &>(*data.parsers[component]); // NOLINT
455  values[component] = parser.Eval();
456  }
457  } // try
458  catch (mu::ParserError &e)
459  {
460  std::cerr << "Message: <" << e.GetMsg() << ">\n";
461  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
462  std::cerr << "Token: <" << e.GetToken() << ">\n";
463  std::cerr << "Position: <" << e.GetPos() << ">\n";
464  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
465  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
466  } // catch
467 #else
468  (void)p;
469  (void)time;
470  (void)values;
472 #endif
473  }
474 
475 // explicit instantiations
476 #include "mu_parser_internal.inst"
477 
478  } // namespace FunctionParser
479 } // namespace internal
480 
Definition: point.h:111
Number do_value(const Point< dim > &p, const double time, unsigned int component) const
void do_all_values(const Point< dim > &p, const double time, ArrayView< Number > &values) const
virtual void initialize(const std::string &vars, const std::vector< std::string > &expressions, const std::map< std::string, double > &constants, const bool time_dependent=false)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:457
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:458
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcParseError(int arg1, std::string arg2)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
static ::ExceptionBase & ExcNeedsFunctionparser()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
Expression ceil(const Expression &x)
Expression floor(const Expression &x)
Expression erfc(const Expression &x)
Expression erf(const Expression &x)
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::vector< std::string > split_string_list(const std::string &s, const std::string &delimiter=",")
Definition: utilities.cc:704
double mu_if(double condition, double thenvalue, double elsevalue)
double mu_erf(double value)
double mu_sec(double value)
double mu_floor(double value)
double mu_csc(double value)
double mu_pow(double a, double b)
double mu_log(double value)
double mu_rand_seed(double seed)
double mu_ceil(double value)
std::vector< std::string > get_function_names()
double mu_cot(double value)
double mu_int(double value)
double mu_or(double left, double right)
double mu_and(double left, double right)
double mu_erfc(double value)
std::vector< std::unique_ptr< muParserBase > > parsers