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