Reference documentation for deal.II version Git 8596a7cd07 2020-12-04 07:30:43 +0100
\(\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\}}\)
function_parser.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2020 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 
16 
19 #include <deal.II/base/patterns.h>
21 #include <deal.II/base/utilities.h>
22 
23 #include <deal.II/lac/vector.h>
24 
25 #include <cmath>
26 #include <map>
27 
28 #ifdef DEAL_II_WITH_MUPARSER
29 # include <muParser.h>
30 #endif
31 
33 
34 
35 template <int dim>
36 const std::vector<std::string> &
38 {
39  return expressions;
40 }
41 
42 
43 
44 template <int dim>
45 FunctionParser<dim>::FunctionParser(const unsigned int n_components,
46  const double initial_time,
47  const double h)
48  : AutoDerivativeFunction<dim>(h, n_components, initial_time)
49  , initialized(false)
50  , n_vars(0)
51 {}
52 
53 
54 template <int dim>
55 FunctionParser<dim>::FunctionParser(const std::string &expression,
56  const std::string &constants,
57  const std::string &variable_names,
58  const double h)
60  h,
61  Utilities::split_string_list(expression, ';').size())
62  , initialized(false)
63  , n_vars(0)
64 {
66  constants,
67  std::make_unique<Patterns::Map>(Patterns::Anything(),
69  0,
71  ",",
72  "="));
73  initialize(variable_names,
74  expression,
75  constants_map,
76  Utilities::split_string_list(variable_names, ",").size() ==
77  dim + 1);
78 }
79 
80 
81 
82 // Note: we explicitly define the destructor here (instead of silently using
83 // the default destructor by declaring nothing in the header) since we do not
84 // expect muParser.h to be available in user projects: i.e., the destructor
85 // must be defined in the source file.
86 template <int dim>
88 
89 
90 #ifdef DEAL_II_WITH_MUPARSER
91 
92 template <int dim>
93 void
94 FunctionParser<dim>::initialize(const std::string & variables,
95  const std::vector<std::string> &expressions,
96  const std::map<std::string, double> &constants,
97  const bool time_dependent)
98 {
99  this->fp.clear(); // this will reset all thread-local objects
100 
101  this->constants = constants;
102  this->var_names = Utilities::split_string_list(variables, ',');
103  this->expressions = expressions;
104  AssertThrow(((time_dependent) ? dim + 1 : dim) == var_names.size(),
105  ExcMessage("Wrong number of variables"));
106 
107  // We check that the number of components of this function matches the
108  // number of components passed in as a vector of strings.
109  AssertThrow(this->n_components == expressions.size(),
110  ExcInvalidExpressionSize(this->n_components, expressions.size()));
111 
112  // Now we define how many variables we expect to read in. We distinguish
113  // between two cases: Time dependent problems, and not time dependent
114  // problems. In the first case the number of variables is given by the
115  // dimension plus one. In the other case, the number of variables is equal
116  // to the dimension. Once we parsed the variables string, if none of this is
117  // the case, then an exception is thrown.
118  if (time_dependent)
119  n_vars = dim + 1;
120  else
121  n_vars = dim;
122 
123  // create a parser object for the current thread we can then query in
124  // value() and vector_value(). this is not strictly necessary because a user
125  // may never call these functions on the current thread, but it gets us
126  // error messages about wrong formulas right away
127  init_muparser();
128 
129  // finally set the initialization bit
130  initialized = true;
131 }
132 
133 
134 
135 template <int dim>
136 void
138 {
139  // check that we have not already initialized the parser on the
140  // current thread, i.e., that the current function is only called
141  // once per thread
142  Assert(fp.get().size() == 0, ExcInternalError());
143 
144  // initialize the objects for the current thread (fp.get() and
145  // vars.get())
146  fp.get().reserve(this->n_components);
147  vars.get().resize(var_names.size());
148  for (unsigned int component = 0; component < this->n_components; ++component)
149  {
150  fp.get().emplace_back(new mu::Parser());
151 
152  for (const auto &constant : constants)
153  {
154  fp.get()[component]->DefineConst(constant.first, constant.second);
155  }
156 
157  for (unsigned int iv = 0; iv < var_names.size(); ++iv)
158  fp.get()[component]->DefineVar(var_names[iv], &vars.get()[iv]);
159 
160  // define some compatibility functions:
161  fp.get()[component]->DefineFun("if",
163  true);
164  fp.get()[component]->DefineOprt("|", internal::FunctionParser::mu_or, 1);
165  fp.get()[component]->DefineOprt("&", internal::FunctionParser::mu_and, 2);
166  fp.get()[component]->DefineFun("int",
168  true);
169  fp.get()[component]->DefineFun("ceil",
171  true);
172  fp.get()[component]->DefineFun("cot",
174  true);
175  fp.get()[component]->DefineFun("csc",
177  true);
178  fp.get()[component]->DefineFun("floor",
180  true);
181  fp.get()[component]->DefineFun("sec",
183  true);
184  fp.get()[component]->DefineFun("log",
186  true);
187  fp.get()[component]->DefineFun("pow",
189  true);
190  fp.get()[component]->DefineFun("erfc",
192  true);
193  fp.get()[component]->DefineFun("rand_seed",
195  true);
196  fp.get()[component]->DefineFun("rand",
198  true);
199 
200  try
201  {
202  // muparser expects that functions have no
203  // space between the name of the function and the opening
204  // parenthesis. this is awkward because it is not backward
205  // compatible to the library we used to use before muparser
206  // (the fparser library) but also makes no real sense.
207  // consequently, in the expressions we set, remove any space
208  // we may find after function names
209  std::string transformed_expression = expressions[component];
210 
211  for (const auto &current_function_name :
213  {
214  const unsigned int function_name_length =
215  current_function_name.size();
216 
217  std::string::size_type pos = 0;
218  while (true)
219  {
220  // try to find any occurrences of the function name
221  pos = transformed_expression.find(current_function_name, pos);
222  if (pos == std::string::npos)
223  break;
224 
225  // replace whitespace until there no longer is any
226  while ((pos + function_name_length <
227  transformed_expression.size()) &&
228  ((transformed_expression[pos + function_name_length] ==
229  ' ') ||
230  (transformed_expression[pos + function_name_length] ==
231  '\t')))
232  transformed_expression.erase(
233  transformed_expression.begin() + pos +
234  function_name_length);
235 
236  // move the current search position by the size of the
237  // actual function name
238  pos += function_name_length;
239  }
240  }
241 
242  // now use the transformed expression
243  fp.get()[component]->SetExpr(transformed_expression);
244  }
245  catch (mu::ParserError &e)
246  {
247  std::cerr << "Message: <" << e.GetMsg() << ">\n";
248  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
249  std::cerr << "Token: <" << e.GetToken() << ">\n";
250  std::cerr << "Position: <" << e.GetPos() << ">\n";
251  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
252  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
253  }
254  }
255 }
256 
257 
258 
259 template <int dim>
260 void
262  const std::string & expression,
263  const std::map<std::string, double> &constants,
264  const bool time_dependent)
265 {
266  initialize(vars,
267  Utilities::split_string_list(expression, ';'),
268  constants,
269  time_dependent);
270 }
271 
272 
273 
274 template <int dim>
275 double
277  const unsigned int component) const
278 {
280  AssertIndexRange(component, this->n_components);
281 
282  // initialize the parser if that hasn't happened yet on the current thread
283  if (fp.get().size() == 0)
284  init_muparser();
285 
286  for (unsigned int i = 0; i < dim; ++i)
287  vars.get()[i] = p(i);
288  if (dim != n_vars)
289  vars.get()[dim] = this->get_time();
290 
291  try
292  {
293  return fp.get()[component]->Eval();
294  }
295  catch (mu::ParserError &e)
296  {
297  std::cerr << "Message: <" << e.GetMsg() << ">\n";
298  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
299  std::cerr << "Token: <" << e.GetToken() << ">\n";
300  std::cerr << "Position: <" << e.GetPos() << ">\n";
301  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
302  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
303  return 0.0;
304  }
305 }
306 
307 
308 
309 template <int dim>
310 void
312  Vector<double> & values) const
313 {
315  Assert(values.size() == this->n_components,
316  ExcDimensionMismatch(values.size(), this->n_components));
317 
318 
319  // initialize the parser if that hasn't happened yet on the current thread
320  if (fp.get().size() == 0)
321  init_muparser();
322 
323  for (unsigned int i = 0; i < dim; ++i)
324  vars.get()[i] = p(i);
325  if (dim != n_vars)
326  vars.get()[dim] = this->get_time();
327 
328  for (unsigned int component = 0; component < this->n_components; ++component)
329  values(component) = fp.get()[component]->Eval();
330 }
331 
332 #else
333 
334 
335 template <int dim>
336 void
337 FunctionParser<dim>::initialize(const std::string &,
338  const std::vector<std::string> &,
339  const std::map<std::string, double> &,
340  const bool)
341 {
343 }
344 
345 template <int dim>
346 void
347 FunctionParser<dim>::initialize(const std::string &,
348  const std::string &,
349  const std::map<std::string, double> &,
350  const bool)
351 {
353 }
354 
355 
356 
357 template <int dim>
358 double
359 FunctionParser<dim>::value(const Point<dim> &, unsigned int) const
360 {
362  return 0.;
363 }
364 
365 
366 template <int dim>
367 void
368 FunctionParser<dim>::vector_value(const Point<dim> &, Vector<double> &) const
369 {
371 }
372 
373 
374 #endif
375 
376 // Explicit Instantiations.
377 
378 template class FunctionParser<1>;
379 template class FunctionParser<2>;
380 template class FunctionParser<3>;
381 
std::vector< std::string > split_string_list(const std::string &s, const std::string &delimiter=",")
Definition: utilities.cc:705
static ::ExceptionBase & ExcNeedsFunctionparser()
double mu_ceil(double value)
FunctionParser(const unsigned int n_components=1, const double initial_time=0.0, const double h=1e-8)
const unsigned int n_components
Definition: function.h:164
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
double mu_int(double value)
static T to_value(const std::string &s, const std::unique_ptr< Patterns::PatternBase > &p=Convert< T >::to_pattern())=delete
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
static const unsigned int max_int_value
Definition: patterns.h:592
unsigned int n_vars
static ::ExceptionBase & ExcParseError(int arg1, std::string arg2)
double mu_floor(double value)
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
double mu_rand_seed(double seed)
double mu_csc(double value)
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const override
double mu_and(double left, double right)
Threads::ThreadLocalStorage< std::vector< std::unique_ptr< mu::Parser > > > fp
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double mu_or(double left, double right)
Threads::ThreadLocalStorage< std::vector< double > > vars
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
void init_muparser() const
static ::ExceptionBase & ExcInvalidExpressionSize(int arg1, int arg2)
double mu_if(double condition, double thenvalue, double elsevalue)
double mu_log(double value)
double mu_pow(double a, double b)
double mu_sec(double value)
double mu_erfc(double value)
Definition: cuda.h:32
std::vector< std::string > var_names
const std::vector< std::string > & get_expressions() const
void initialize(const std::string &vars, const std::vector< std::string > &expressions, const ConstMap &constants, const bool time_dependent=false)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
Number get_time() const
std::vector< std::string > function_names
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
double mu_cot(double value)
std::vector< std::string > expressions
virtual ~FunctionParser() override
std::map< std::string, double > constants
static ::ExceptionBase & ExcInternalError()