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