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