Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
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
34namespace 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
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
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>
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,
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()),
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()),
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: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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNeedsFunctionparser()
#define Assert(cond, exc)
static ::ExceptionBase & ExcParseError(int arg1, std::string arg2)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
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)
::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)
std::vector< std::unique_ptr< muParserBase > > parsers