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\}}\)
tensor_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 
18 #include <deal.II/base/patterns.h>
19 #include <deal.II/base/tensor.h>
23 #include <deal.II/base/utilities.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 rank, int dim, typename Number>
36 const std::vector<std::string> &
38 {
39  return expressions;
40 }
41 
42 
43 
44 template <int rank, int dim, typename Number>
46  const double initial_time)
47  : TensorFunction<rank, dim, Number>(initial_time)
48  , initialized(false)
49  , n_vars(0)
50  , n_components(Utilities::pow(dim, rank))
51 {}
52 
53 
54 template <int rank, int dim, typename Number>
56  const std::string &expression,
57  const std::string &constants,
58  const std::string &variable_names)
59  : TensorFunction<rank, dim, Number>()
60  , initialized(false)
61  , n_vars(0)
62  , n_components(Utilities::pow(dim, rank))
63 {
65  constants,
66  std::make_unique<Patterns::Map>(Patterns::Anything(),
68  0,
70  ",",
71  "="));
72  initialize(variable_names,
73  expression,
74  constants_map,
75  Utilities::split_string_list(variable_names, ",").size() ==
76  dim + 1);
77 }
78 
79 
80 // We deliberately delay the definition of the default destructor
81 // so that we don't need to include the definition of mu::Parser
82 // in the header file.
83 template <int rank, int dim, typename Number>
85 
86 
87 #ifdef DEAL_II_WITH_MUPARSER
88 
89 template <int rank, int dim, typename Number>
90 void
92  const std::string & variables,
93  const std::vector<std::string> & expressions,
94  const std::map<std::string, double> &constants,
95  const bool time_dependent)
96 {
97  this->tfp.clear(); // this will reset all thread-local objects
98 
99  this->constants = constants;
100  this->var_names = Utilities::split_string_list(variables, ',');
101  this->expressions = expressions;
102  AssertThrow(((time_dependent) ? dim + 1 : dim) == var_names.size(),
103  ExcMessage("Wrong number of variables"));
104 
105  // We check that the number of
106  // components of this function
107  // matches the number of components
108  // passed in as a vector of
109  // strings.
110  AssertThrow(this->n_components == expressions.size(),
111  ExcInvalidExpressionSize(this->n_components, expressions.size()));
112 
113  // Now we define how many variables
114  // we expect to read in. We
115  // distinguish between two cases:
116  // Time dependent problems, and not
117  // time dependent problems. In the
118  // first case the number of
119  // variables is given by the
120  // dimension plus one. In the other
121  // case, the number of variables is
122  // equal to the dimension. Once we
123  // parsed the variables string, if
124  // none of this is the case, then
125  // an exception is thrown.
126  if (time_dependent)
127  n_vars = dim + 1;
128  else
129  n_vars = dim;
130 
131  // create a parser object for the current thread we can then query
132  // in value() and vector_value(). this is not strictly necessary
133  // because a user may never call these functions on the current
134  // thread, but it gets us error messages about wrong formulas right
135  // away
136  init_muparser();
137 
138  // finally set the initialization bit
139  initialized = true;
140 }
141 
142 
143 
144 template <int rank, int dim, typename Number>
145 void
147 {
148  // check that we have not already initialized the parser on the
149  // current thread, i.e., that the current function is only called
150  // once per thread
151  Assert(tfp.get().size() == 0, ExcInternalError());
152 
153  // initialize the objects for the current thread (tfp.get() and
154  // vars.get())
155  tfp.get().reserve(this->n_components);
156  vars.get().resize(var_names.size());
157  for (unsigned int component = 0; component < this->n_components; ++component)
158  {
159  tfp.get().emplace_back(new mu::Parser());
160 
161  for (const auto &constant : constants)
162  {
163  tfp.get()[component]->DefineConst(constant.first, constant.second);
164  }
165 
166  for (unsigned int iv = 0; iv < var_names.size(); ++iv)
167  tfp.get()[component]->DefineVar(var_names[iv], &vars.get()[iv]);
168 
169  // define some compatibility functions:
170  tfp.get()[component]->DefineFun("if",
172  true);
173  tfp.get()[component]->DefineOprt("|", internal::FunctionParser::mu_or, 1);
174  tfp.get()[component]->DefineOprt("&",
176  2);
177  tfp.get()[component]->DefineFun("int",
179  true);
180  tfp.get()[component]->DefineFun("ceil",
182  true);
183  tfp.get()[component]->DefineFun("cot",
185  true);
186  tfp.get()[component]->DefineFun("csc",
188  true);
189  tfp.get()[component]->DefineFun("floor",
191  true);
192  tfp.get()[component]->DefineFun("sec",
194  true);
195  tfp.get()[component]->DefineFun("log",
197  true);
198  tfp.get()[component]->DefineFun("pow",
200  true);
201  tfp.get()[component]->DefineFun("erfc",
203  true);
204  tfp.get()[component]->DefineFun("rand_seed",
206  true);
207  tfp.get()[component]->DefineFun("rand",
209  true);
210 
211  try
212  {
213  // muparser expects that functions have no
214  // space between the name of the function and the opening
215  // parenthesis. this is awkward because it is not backward
216  // compatible to the library we used to use before muparser
217  // (the tfparser library) but also makes no real sense.
218  // consequently, in the expressions we set, remove any space
219  // we may find after function names
220  std::string transformed_expression = expressions[component];
221 
222  for (const auto &current_function_name :
224  {
225  const unsigned int function_name_length =
226  current_function_name.size();
227 
228  std::string::size_type pos = 0;
229  while (true)
230  {
231  // try to find any occurrences of the function name
232  pos = transformed_expression.find(current_function_name, pos);
233  if (pos == std::string::npos)
234  break;
235 
236  // replace whitespace until there no longer is any
237  while ((pos + function_name_length <
238  transformed_expression.size()) &&
239  ((transformed_expression[pos + function_name_length] ==
240  ' ') ||
241  (transformed_expression[pos + function_name_length] ==
242  '\t')))
243  transformed_expression.erase(
244  transformed_expression.begin() + pos +
245  function_name_length);
246 
247  // move the current search position by the size of the
248  // actual function name
249  pos += function_name_length;
250  }
251  }
252 
253  // now use the transformed expression
254  tfp.get()[component]->SetExpr(transformed_expression);
255  }
256  catch (mu::ParserError &e)
257  {
258  std::cerr << "Message: <" << e.GetMsg() << ">\n";
259  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
260  std::cerr << "Token: <" << e.GetToken() << ">\n";
261  std::cerr << "Position: <" << e.GetPos() << ">\n";
262  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
263  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
264  }
265  }
266 }
267 
268 
269 
270 template <int rank, int dim, typename Number>
271 void
273  const std::string & vars,
274  const std::string & expression,
275  const std::map<std::string, double> &constants,
276  const bool time_dependent)
277 {
278  initialize(vars,
279  Utilities::split_string_list(expression, ';'),
280  constants,
281  time_dependent);
282 }
283 
284 
285 
286 template <int rank, int dim, typename Number>
289 {
291 
292  // initialize the parser if that hasn't happened yet on the current thread
293  if (tfp.get().size() == 0)
294  init_muparser();
295 
296  for (unsigned int i = 0; i < dim; ++i)
297  vars.get()[i] = p(i);
298  if (dim != n_vars)
299  vars.get()[dim] = this->get_time();
300 
301  // initialize tensor with zeros
303 
304  try
305  {
306  unsigned int component = 0;
307  for (Number *value_ptr = value.begin_raw(); value_ptr != value.end_raw();
308  ++value_ptr)
309  {
310  *value_ptr = tfp.get()[component]->Eval();
311  ++component;
312  } // for
313  } // try
314  catch (mu::ParserError &e)
315  {
316  std::cerr << "Message: <" << e.GetMsg() << ">\n";
317  std::cerr << "Formula: <" << e.GetExpr() << ">\n";
318  std::cerr << "Token: <" << e.GetToken() << ">\n";
319  std::cerr << "Position: <" << e.GetPos() << ">\n";
320  std::cerr << "Errc: <" << e.GetCode() << ">" << std::endl;
321  AssertThrow(false, ExcParseError(e.GetCode(), e.GetMsg()));
322 
323  return Tensor<rank, dim, Number>();
324  } // catch
325 
326  return value;
327 }
328 
329 
330 
331 template <int rank, int dim, typename Number>
332 void
334  const std::vector<Point<dim>> & p,
335  std::vector<Tensor<rank, dim, Number>> &values) const
336 {
337  Assert(p.size() == values.size(),
338  ExcDimensionMismatch(p.size(), values.size()));
339 
340  for (unsigned int i = 0; i < p.size(); ++i)
341  {
342  values[i] = value(p[i]);
343  }
344 }
345 
346 #else
347 
348 
349 template <int rank, int dim, typename Number>
350 void
352  const std::string &,
353  const std::vector<std::string> &,
354  const std::map<std::string, double> &,
355  const bool)
356 {
358 }
359 
360 template <int rank, int dim, typename Number>
361 void
363  const std::string &,
364  const std::string &,
365  const std::map<std::string, double> &,
366  const bool)
367 {
369 }
370 
371 
372 
373 template <int rank, int dim, typename Number>
376 {
378  return Tensor<rank, dim, Number>();
379 }
380 
381 
382 
383 template <int rank, int dim, typename Number>
384 void
386  const std::vector<Point<dim>> &,
387  std::vector<Tensor<rank, dim, Number>> &) const
388 {
390 }
391 
392 
393 #endif
394 
395 // explicit instantiations
396 #include "tensor_function_parser.inst"
397 
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)
static ::ExceptionBase & ExcParseError(int arg1, std::string arg2)
Number * begin_raw()
std::map< std::string, double > constants
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
static const unsigned int max_int_value
Definition: patterns.h:592
Threads::ThreadLocalStorage< std::vector< std::unique_ptr< mu::Parser > > > tfp
double mu_floor(double value)
virtual ~TensorFunctionParser() override
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
double mu_rand_seed(double seed)
double mu_csc(double value)
virtual Tensor< rank, dim, Number > value(const Point< dim > &p) const override
Number * end_raw()
static ::ExceptionBase & ExcInvalidExpressionSize(int arg1, int arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
double mu_and(double left, double right)
#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)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &x, const Number p)
double mu_if(double condition, double thenvalue, double elsevalue)
double mu_log(double value)
void initialize(const std::string &vars, const std::vector< std::string > &expressions, const ConstMap &constants, const bool time_dependent=false)
double mu_pow(double a, double b)
TensorFunctionParser(const double initial_time=0.0)
virtual void value_list(const std::vector< Point< dim >> &p, std::vector< Tensor< rank, dim, Number >> &values) const override
double mu_sec(double value)
double mu_erfc(double value)
Definition: cuda.h:32
std::vector< std::string > expressions
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
Threads::ThreadLocalStorage< std::vector< double > > vars
numbers::NumberTraits< Number >::real_type get_time() const
std::vector< std::string > function_names
std::vector< std::string > var_names
double mu_cot(double value)
const std::vector< std::string > & get_expressions() const
static ::ExceptionBase & ExcInternalError()