deal.II version GIT relicensing-1936-ge55124254b 2024-10-04 21:10:00+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 return uniform_distribution(
171 rng_map.try_emplace(seed, std::mt19937(static_cast<unsigned int>(seed)))
172 .first->second);
173 }
174
175
176 // Returns a random value in the range [0,1]
177 double
179 {
180 static std::mutex rand_mutex;
181 std::lock_guard<std::mutex> lock(rand_mutex);
182 std::uniform_real_distribution<> uniform_distribution(0., 1.);
183 const unsigned int seed = static_cast<unsigned long>(std::time(nullptr));
184 static std::mt19937 rng(seed);
185 return uniform_distribution(rng);
187
188
189
190 std::vector<std::string>
192 {
193 return {// functions predefined by muparser
194 "sin",
195 "cos",
196 "tan",
197 "asin",
198 "acos",
199 "atan",
200 "sinh",
201 "cosh",
202 "tanh",
203 "asinh",
204 "acosh",
205 "atanh",
206 "atan2",
207 "log2",
208 "log10",
209 "log",
210 "ln",
211 "exp",
212 "sqrt",
213 "sign",
214 "rint",
215 "abs",
216 "min",
217 "max",
218 "sum",
219 "avg",
220 // functions we define ourselves above
221 "if",
222 "int",
223 "ceil",
224 "cot",
225 "csc",
226 "floor",
227 "sec",
228 "pow",
229 "erf",
230 "erfc",
231 "rand",
232 "rand_seed"};
233 }
234
235#ifdef DEAL_II_WITH_MUPARSER
239 class Parser : public muParserBase
240 {
241 public:
242 operator mu::Parser &()
243 {
244 return parser;
245 }
246
247 operator const mu::Parser &() const
248 {
249 return parser;
250 }
251
252 protected:
253 mu::Parser parser;
254 };
255#endif
256
257
258
259 template <int dim, typename Number>
261 : initialized(false)
262 , n_vars(0)
263 {}
264
265
266
267 template <int dim, typename Number>
268 void
270 const std::string &variables,
271 const std::vector<std::string> &expressions,
272 const std::map<std::string, double> &constants,
273 const bool time_dependent)
274 {
275 this->parser_data.clear(); // this will reset all thread-local objects
276
277 this->constants = constants;
278 this->var_names = Utilities::split_string_list(variables, ',');
279 this->expressions = expressions;
280 AssertThrow(((time_dependent) ? dim + 1 : dim) == this->var_names.size(),
281 ExcMessage("Wrong number of variables"));
282
283 // Now we define how many variables we expect to read in. We distinguish
284 // between two cases: Time dependent problems, and not time dependent
285 // problems. In the first case the number of variables is given by the
286 // dimension plus one. In the other case, the number of variables is equal
287 // to the dimension. Once we parsed the variables string, if none of this
288 // is the case, then an exception is thrown.
289 if (time_dependent)
290 this->n_vars = dim + 1;
291 else
292 this->n_vars = dim;
293
294 // create a parser object for the current thread we can then query in
295 // value() and vector_value(). this is not strictly necessary because a
296 // user may never call these functions on the current thread, but it gets
297 // us error messages about wrong formulas right away
298 this->init_muparser();
299 this->initialized = true;
300 }
301
302
303
304 template <int dim, typename Number>
305 void
307 {
308#ifdef DEAL_II_WITH_MUPARSER
309 // check that we have not already initialized the parser on the
310 // current thread, i.e., that the current function is only called
311 // once per thread
312 ParserData &data = this->parser_data.get();
313 Assert(data.parsers.empty() && data.vars.empty(), ExcInternalError());
314 const unsigned int n_components = expressions.size();
315
316 // initialize the objects for the current thread
317 data.parsers.reserve(n_components);
318 data.vars.resize(this->var_names.size());
319 for (unsigned int component = 0; component < n_components; ++component)
320 {
321 data.parsers.emplace_back(std::make_unique<Parser>());
322 mu::Parser &parser = dynamic_cast<Parser &>(*data.parsers.back());
323
324 for (const auto &constant : this->constants)
325 parser.DefineConst(constant.first, constant.second);
326
327 for (unsigned int iv = 0; iv < this->var_names.size(); ++iv)
328 parser.DefineVar(this->var_names[iv], &data.vars[iv]);
329
330 // define some compatibility functions:
331 parser.DefineFun("if", mu_if, true);
332 parser.DefineOprt("|", mu_or, 1);
333 parser.DefineOprt("&", mu_and, 2);
334 parser.DefineFun("int", mu_int, true);
335 parser.DefineFun("ceil", mu_ceil, true);
336 parser.DefineFun("cot", mu_cot, true);
337 parser.DefineFun("csc", mu_csc, true);
338 parser.DefineFun("floor", mu_floor, true);
339 parser.DefineFun("sec", mu_sec, true);
340 parser.DefineFun("log", mu_log, true);
341 parser.DefineFun("pow", mu_pow, true);
342 parser.DefineFun("erfc", mu_erfc, true);
343 // Disable optimizations (by passing false) that assume the functions
344 // will always return the same value:
345 parser.DefineFun("rand_seed", mu_rand_seed, false);
346 parser.DefineFun("rand", mu_rand, false);
347
348 try
349 {
350 // muparser expects that functions have no
351 // space between the name of the function and the opening
352 // parenthesis. this is awkward because it is not backward
353 // compatible to the library we used to use before muparser
354 // (the fparser library) but also makes no real sense.
355 // consequently, in the expressions we set, remove any space
356 // we may find after function names
357 std::string transformed_expression = this->expressions[component];
358
359 for (const auto &current_function_name : get_function_names())
360 {
361 const unsigned int function_name_length =
362 current_function_name.size();
363
364 std::string::size_type pos = 0;
365 while (true)
366 {
367 // try to find any occurrences of the function name
368 pos =
369 transformed_expression.find(current_function_name, pos);
370 if (pos == std::string::npos)
371 break;
372
373 // replace whitespace until there no longer is any
374 while (
375 (pos + function_name_length <
376 transformed_expression.size()) &&
377 ((transformed_expression[pos + function_name_length] ==
378 ' ') ||
379 (transformed_expression[pos + function_name_length] ==
380 '\t')))
381 transformed_expression.erase(
382 transformed_expression.begin() + pos +
383 function_name_length);
384
385 // move the current search position by the size of the
386 // actual function name
387 pos += function_name_length;
388 }
389 }
390
391 // now use the transformed expression
392 parser.SetExpr(transformed_expression);
393 }
394 catch (mu::ParserError &e)
395 {
396 std::cerr << "Message: <" << e.GetMsg() << ">\n";
397 std::cerr << "Formula: <" << e.GetExpr() << ">\n";
398 std::cerr << "Token: <" << e.GetToken() << ">\n";
399 std::cerr << "Position: <" << e.GetPos() << ">\n";
400 std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
401 AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
402 }
403 }
404#else
406#endif
407 }
408
409 template <int dim, typename Number>
410 Number
412 const double time,
413 unsigned int component) const
414 {
415#ifdef DEAL_II_WITH_MUPARSER
416 Assert(this->initialized == true, ExcNotInitialized());
417
418 // initialize the parser if that hasn't happened yet on the current
419 // thread
420 internal::FunctionParser::ParserData &data = this->parser_data.get();
421 if (data.vars.empty())
422 init_muparser();
423
424 for (unsigned int i = 0; i < dim; ++i)
425 data.vars[i] = p[i];
426 if (dim != this->n_vars)
427 data.vars[dim] = time;
428
429 try
430 {
431 Assert(dynamic_cast<Parser *>(data.parsers[component].get()),
433 // NOLINTNEXTLINE don't warn about using static_cast once we check
434 mu::Parser &parser = static_cast<Parser &>(*data.parsers[component]);
435 return parser.Eval();
436 } // try
437 catch (mu::ParserError &e)
438 {
439 std::cerr << "Message: <" << e.GetMsg() << ">\n";
440 std::cerr << "Formula: <" << e.GetExpr() << ">\n";
441 std::cerr << "Token: <" << e.GetToken() << ">\n";
442 std::cerr << "Position: <" << e.GetPos() << ">\n";
443 std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
444 AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
445 } // catch
446
447#else
448 (void)p;
449 (void)time;
450 (void)component;
452#endif
453 return std::numeric_limits<double>::signaling_NaN();
454 }
455
456 template <int dim, typename Number>
457 void
459 const Point<dim> &p,
460 const double time,
461 ArrayView<Number> &values) const
462 {
463#ifdef DEAL_II_WITH_MUPARSER
464 Assert(this->initialized == true, ExcNotInitialized());
465
466 // initialize the parser if that hasn't happened yet on the current
467 // thread
468 internal::FunctionParser::ParserData &data = this->parser_data.get();
469 if (data.vars.empty())
470 init_muparser();
471
472 for (unsigned int i = 0; i < dim; ++i)
473 data.vars[i] = p[i];
474 if (dim != this->n_vars)
475 data.vars[dim] = time;
476
477 AssertDimension(values.size(), data.parsers.size());
478 try
479 {
480 for (unsigned int component = 0; component < data.parsers.size();
481 ++component)
482 {
483 Assert(dynamic_cast<Parser *>(data.parsers[component].get()),
485 mu::Parser &parser =
486 // We just checked that the pointer is valid so suppress the
487 // clang-tidy check
488 static_cast<Parser &>(*data.parsers[component]); // NOLINT
489 values[component] = parser.Eval();
490 }
491 } // try
492 catch (mu::ParserError &e)
493 {
494 std::cerr << "Message: <" << e.GetMsg() << ">\n";
495 std::cerr << "Formula: <" << e.GetExpr() << ">\n";
496 std::cerr << "Token: <" << e.GetToken() << ">\n";
497 std::cerr << "Position: <" << e.GetPos() << ">\n";
498 std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
499 AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
500 } // catch
501#else
502 (void)p;
503 (void)time;
504 (void)values;
506#endif
507 }
508
509// explicit instantiations
510#include "mu_parser_internal.inst"
511
512 } // namespace FunctionParser
513} // namespace internal
514
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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