Reference documentation for deal.II version GIT 67d6815257 2023-03-20 17:20:01+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  // Disable optimizations (by passing false) that assume the functions
312  // will always return the same value:
313  parser.DefineFun("rand_seed", mu_rand_seed, false);
314  parser.DefineFun("rand", mu_rand, false);
315 
316  try
317  {
318  // muparser expects that functions have no
319  // space between the name of the function and the opening
320  // parenthesis. this is awkward because it is not backward
321  // compatible to the library we used to use before muparser
322  // (the fparser library) but also makes no real sense.
323  // consequently, in the expressions we set, remove any space
324  // we may find after function names
325  std::string transformed_expression = this->expressions[component];
326 
327  for (const auto &current_function_name : get_function_names())
328  {
329  const unsigned int function_name_length =
330  current_function_name.size();
331 
332  std::string::size_type pos = 0;
333  while (true)
334  {
335  // try to find any occurrences of the function name
336  pos =
337  transformed_expression.find(current_function_name, pos);
338  if (pos == std::string::npos)
339  break;
340 
341  // replace whitespace until there no longer is any
342  while (
343  (pos + function_name_length <
344  transformed_expression.size()) &&
345  ((transformed_expression[pos + function_name_length] ==
346  ' ') ||
347  (transformed_expression[pos + function_name_length] ==
348  '\t')))
349  transformed_expression.erase(
350  transformed_expression.begin() + pos +
351  function_name_length);
352 
353  // move the current search position by the size of the
354  // actual function name
355  pos += function_name_length;
356  }
357  }
358 
359  // now use the transformed expression
360  parser.SetExpr(transformed_expression);
361  }
362  catch (mu::ParserError &e)
363  {
364  std::cerr << "Message: <" << e.GetMsg() << ">\n";
365  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
366  std::cerr << "Token: <" << e.GetToken() << ">\n";
367  std::cerr << "Position: <" << e.GetPos() << ">\n";
368  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
369  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
370  }
371  }
372 #else
374 #endif
375  }
376 
377  template <int dim, typename Number>
378  Number
380  const double time,
381  unsigned int component) const
382  {
383 #ifdef DEAL_II_WITH_MUPARSER
384  Assert(this->initialized == true, ExcNotInitialized());
385 
386  // initialize the parser if that hasn't happened yet on the current
387  // thread
388  internal::FunctionParser::ParserData &data = this->parser_data.get();
389  if (data.vars.size() == 0)
390  init_muparser();
391 
392  for (unsigned int i = 0; i < dim; ++i)
393  data.vars[i] = p(i);
394  if (dim != this->n_vars)
395  data.vars[dim] = time;
396 
397  try
398  {
399  Assert(dynamic_cast<Parser *>(data.parsers[component].get()),
400  ExcInternalError());
401  // NOLINTNEXTLINE don't warn about using static_cast once we check
402  mu::Parser &parser = static_cast<Parser &>(*data.parsers[component]);
403  return parser.Eval();
404  } // try
405  catch (mu::ParserError &e)
406  {
407  std::cerr << "Message: <" << e.GetMsg() << ">\n";
408  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
409  std::cerr << "Token: <" << e.GetToken() << ">\n";
410  std::cerr << "Position: <" << e.GetPos() << ">\n";
411  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
412  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
413  } // catch
414 
415 #else
416  (void)p;
417  (void)time;
418  (void)component;
420 #endif
421  return std::numeric_limits<double>::signaling_NaN();
422  }
423 
424  template <int dim, typename Number>
425  void
427  const Point<dim> & p,
428  const double time,
429  ArrayView<Number> &values) const
430  {
431 #ifdef DEAL_II_WITH_MUPARSER
432  Assert(this->initialized == true, ExcNotInitialized());
433 
434  // initialize the parser if that hasn't happened yet on the current
435  // thread
436  internal::FunctionParser::ParserData &data = this->parser_data.get();
437  if (data.vars.size() == 0)
438  init_muparser();
439 
440  for (unsigned int i = 0; i < dim; ++i)
441  data.vars[i] = p(i);
442  if (dim != this->n_vars)
443  data.vars[dim] = time;
444 
445  AssertDimension(values.size(), data.parsers.size());
446  try
447  {
448  for (unsigned int component = 0; component < data.parsers.size();
449  ++component)
450  {
451  Assert(dynamic_cast<Parser *>(data.parsers[component].get()),
452  ExcInternalError());
453  mu::Parser &parser =
454  // We just checked that the pointer is valid so suppress the
455  // clang-tidy check
456  static_cast<Parser &>(*data.parsers[component]); // NOLINT
457  values[component] = parser.Eval();
458  }
459  } // try
460  catch (mu::ParserError &e)
461  {
462  std::cerr << "Message: <" << e.GetMsg() << ">\n";
463  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
464  std::cerr << "Token: <" << e.GetToken() << ">\n";
465  std::cerr << "Position: <" << e.GetPos() << ">\n";
466  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
467  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
468  } // catch
469 #else
470  (void)p;
471  (void)time;
472  (void)values;
474 #endif
475  }
476 
477 // explicit instantiations
478 #include "mu_parser_internal.inst"
479 
480  } // namespace FunctionParser
481 } // namespace internal
482 
Definition: point.h:110
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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcParseError(int arg1, std::string arg2)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcNeedsFunctionparser()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
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