libMesh
parsed_function.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #ifndef LIBMESH_PARSED_FUNCTION_H
19 #define LIBMESH_PARSED_FUNCTION_H
20 
21 #include "libmesh/libmesh_config.h"
22 #include "libmesh/function_base.h"
23 
24 #ifdef LIBMESH_HAVE_FPARSER
25 
26 // Local includes
27 #include "libmesh/dense_vector.h"
28 #include "libmesh/vector_value.h"
29 #include "libmesh/point.h"
30 
31 // FParser includes
32 #include "libmesh/fparser_ad.hh"
33 
34 // C++ includes
35 #include <algorithm> // std::find
36 #include <cmath>
37 #include <cmath>
38 #include <cstddef>
39 #include <iomanip>
40 #include <limits>
41 #include <sstream>
42 #include <string>
43 #include <vector>
44 
45 namespace libMesh
46 {
47 
57 template <typename Output=Number, typename OutputGradient=Gradient>
58 class ParsedFunction : public FunctionBase<Output>
59 {
60 public:
61  explicit
62  ParsedFunction (const std::string & expression,
63  const std::vector<std::string> * additional_vars=libmesh_nullptr,
64  const std::vector<Output> * initial_vals=libmesh_nullptr);
65 
69  void reparse (const std::string & expression);
70 
71  virtual Output operator() (const Point & p,
72  const Real time = 0) libmesh_override;
73 
77  virtual bool has_derivatives() { return _valid_derivatives; }
78 
79  virtual Output dot(const Point & p,
80  const Real time = 0);
81 
82  virtual OutputGradient gradient(const Point & p,
83  const Real time = 0);
84 
85  virtual void operator() (const Point & p,
86  const Real time,
87  DenseVector<Output> & output) libmesh_override;
88 
89  virtual Output component (unsigned int i,
90  const Point & p,
91  Real time) libmesh_override;
92 
93  const std::string & expression() { return _expression; }
94 
98  virtual Output & getVarAddress(const std::string & variable_name);
99 
100  virtual UniquePtr<FunctionBase<Output>> clone() const libmesh_override;
101 
110  Output get_inline_value(const std::string & inline_var_name) const;
111 
122  void set_inline_value(const std::string & inline_var_name,
123  Output newval);
124 
125 protected:
129  void partial_reparse (const std::string & expression);
130 
134  std::size_t find_name (const std::string & varname,
135  const std::string & expr) const;
136 
140  bool expression_is_time_dependent( const std::string & expression ) const;
141 
142 private:
146  void set_spacetime(const Point & p,
147  const Real time = 0);
148 
152  inline Output eval(FunctionParserADBase<Output> & parser,
153  const std::string & libmesh_dbg_var(function_name),
154  unsigned int libmesh_dbg_var(component_idx)) const;
155 
156  std::string _expression;
157  std::vector<std::string> _subexpressions;
158  std::vector<FunctionParserADBase<Output>> parsers;
159  std::vector<Output> _spacetime;
160 
161  // derivative functions
162  std::vector<FunctionParserADBase<Output>> dx_parsers;
163 #if LIBMESH_DIM > 1
164  std::vector<FunctionParserADBase<Output>> dy_parsers;
165 #endif
166 #if LIBMESH_DIM > 2
167  std::vector<FunctionParserADBase<Output>> dz_parsers;
168 #endif
169  std::vector<FunctionParserADBase<Output>> dt_parsers;
171 
172  // Variables/values that can be parsed and handled by the function parser
173  std::string variables;
174  std::vector<std::string> _additional_vars;
175  std::vector<Output> _initial_vals;
176 };
177 
178 
179 /*----------------------- Inline functions ----------------------------------*/
180 
181 template <typename Output, typename OutputGradient>
182 inline
184  const std::vector<std::string> * additional_vars,
185  const std::vector<Output> * initial_vals) :
186  _expression (), // overridden by parse()
187  // Size the spacetime vector to account for space, time, and any additional
188  // variables passed
189  _spacetime (LIBMESH_DIM+1 + (additional_vars ? additional_vars->size() : 0)),
190  _valid_derivatives (true),
191  _additional_vars (additional_vars ? *additional_vars : std::vector<std::string>()),
192  _initial_vals (initial_vals ? *initial_vals : std::vector<Output>())
193 {
194  // time-dependence established in reparse function
195  this->reparse(expression);
196  this->_initialized = true;
197 }
198 
199 
200 template <typename Output, typename OutputGradient>
201 inline
202 void
204 {
205  variables = "x";
206 #if LIBMESH_DIM > 1
207  variables += ",y";
208 #endif
209 #if LIBMESH_DIM > 2
210  variables += ",z";
211 #endif
212  variables += ",t";
213 
214  // If additional vars were passed, append them to the string
215  // that we send to the function parser. Also add them to the
216  // end of our spacetime vector
217  for (std::size_t i=0; i < _additional_vars.size(); ++i)
218  {
219  variables += "," + _additional_vars[i];
220  // Initialize extra variables to the vector passed in or zero
221  // Note: The initial_vals vector can be shorter than the additional_vars vector
222  _spacetime[LIBMESH_DIM+1 + i] =
223  (i < _initial_vals.size()) ? _initial_vals[i] : 0;
224  }
225 
226  this->partial_reparse(expression);
227 
228  this->_is_time_dependent = this->expression_is_time_dependent(expression);
229 }
230 
231 template <typename Output, typename OutputGradient>
232 inline
233 Output
235 {
236  set_spacetime(p, time);
237  return eval(parsers[0], "f", 0);
238 }
239 
240 template <typename Output, typename OutputGradient>
241 inline
242 Output
244 {
245  set_spacetime(p, time);
246  return eval(dt_parsers[0], "df/dt", 0);
247 }
248 
249 template <typename Output, typename OutputGradient>
250 inline
251 OutputGradient
253 {
254  OutputGradient grad;
255  set_spacetime(p, time);
256 
257  grad(0) = eval(dx_parsers[0], "df/dx", 0);
258 #if LIBMESH_DIM > 1
259  grad(1) = eval(dy_parsers[0], "df/dy", 0);
260 #endif
261 #if LIBMESH_DIM > 2
262  grad(2) = eval(dz_parsers[0], "df/dz", 0);
263 #endif
264 
265  return grad;
266 }
267 
268 template <typename Output, typename OutputGradient>
269 inline
270 void
272  (const Point & p,
273  const Real time,
274  DenseVector<Output> & output)
275 {
276  set_spacetime(p, time);
277 
278  unsigned int size = output.size();
279 
280  libmesh_assert_equal_to (size, parsers.size());
281 
282  // The remaining locations in _spacetime are currently fixed at construction
283  // but could potentially be made dynamic
284  for (unsigned int i=0; i != size; ++i)
285  output(i) = eval(parsers[i], "f", i);
286 }
287 
292 template <typename Output, typename OutputGradient>
293 inline
294 Output
296  const Point & p,
297  Real time)
298 {
299  set_spacetime(p, time);
300  libmesh_assert_less (i, parsers.size());
301 
302  // The remaining locations in _spacetime are currently fixed at construction
303  // but could potentially be made dynamic
304  libmesh_assert_less(i, parsers.size());
305  return eval(parsers[i], "f", i);
306 }
307 
311 template <typename Output, typename OutputGradient>
312 inline
313 Output &
314 ParsedFunction<Output,OutputGradient>::getVarAddress (const std::string & variable_name)
315 {
316  const std::vector<std::string>::iterator it =
317  std::find(_additional_vars.begin(), _additional_vars.end(), variable_name);
318 
319  if (it == _additional_vars.end())
320  libmesh_error_msg("ERROR: Requested variable not found in parsed function");
321 
322  // Iterator Arithmetic (How far from the end of the array is our target address?)
323  return _spacetime[_spacetime.size() - (_additional_vars.end() - it)];
324 }
325 
326 
327 template <typename Output, typename OutputGradient>
328 inline
331 {
334 }
335 
336 template <typename Output, typename OutputGradient>
337 inline
338 Output
339 ParsedFunction<Output,OutputGradient>::get_inline_value (const std::string & inline_var_name) const
340 {
341  libmesh_assert_greater (_subexpressions.size(), 0);
342 
343 #ifndef NDEBUG
344  bool found_var_name = false;
345 #endif
346  Output old_var_value(0.);
347 
348  for (std::size_t s=0; s != _subexpressions.size(); ++s)
349  {
350  const std::string & subexpression = _subexpressions[s];
351  const std::size_t varname_i =
352  find_name(inline_var_name, subexpression);
353  if (varname_i == std::string::npos)
354  continue;
355 
356  const std::size_t assignment_i =
357  subexpression.find(":", varname_i+1);
358 
359  libmesh_assert_not_equal_to(assignment_i, std::string::npos);
360 
361  libmesh_assert_equal_to(subexpression[assignment_i+1], '=');
362  for (unsigned int i = varname_i+1; i != assignment_i; ++i)
363  libmesh_assert_equal_to(subexpression[i], ' ');
364 
365  std::size_t end_assignment_i =
366  subexpression.find(";", assignment_i+1);
367 
368  libmesh_assert_not_equal_to(end_assignment_i, std::string::npos);
369 
370  std::string new_subexpression =
371  subexpression.substr(0, end_assignment_i+1) +
372  inline_var_name;
373 
374 #ifdef LIBMESH_HAVE_FPARSER
375  // Parse and evaluate the new subexpression.
376  // Add the same constants as we used originally.
377  FunctionParserADBase<Output> fp;
378  fp.AddConstant("NaN", std::numeric_limits<Real>::quiet_NaN());
379  fp.AddConstant("pi", std::acos(Real(-1)));
380  fp.AddConstant("e", std::exp(Real(1)));
381  if (fp.Parse(new_subexpression, variables) != -1) // -1 for success
382  libmesh_error_msg
383  ("ERROR: FunctionParser is unable to parse modified expression: "
384  << new_subexpression << '\n' << fp.ErrorMsg());
385 
386  Output new_var_value = this->eval(fp, new_subexpression, 0);
387 #ifdef NDEBUG
388  return new_var_value;
389 #else
390  if (found_var_name)
391  {
392  libmesh_assert_equal_to(old_var_value, new_var_value);
393  }
394  else
395  {
396  old_var_value = new_var_value;
397  found_var_name = true;
398  }
399 #endif
400 
401 #else
402  libmesh_error_msg("ERROR: This functionality requires fparser!");
403 #endif
404  }
405 
406  libmesh_assert(found_var_name);
407  return old_var_value;
408 }
409 
410 
411 template <typename Output, typename OutputGradient>
412 inline
413 void
414 ParsedFunction<Output,OutputGradient>::set_inline_value (const std::string & inline_var_name,
415  Output newval)
416 {
417  libmesh_assert_greater (_subexpressions.size(), 0);
418 
419 #ifndef NDEBUG
420  bool found_var_name = false;
421 #endif
422  for (std::size_t s=0; s != _subexpressions.size(); ++s)
423  {
424  const std::string & subexpression = _subexpressions[s];
425  const std::size_t varname_i =
426  find_name(inline_var_name, subexpression);
427  if (varname_i == std::string::npos)
428  continue;
429 
430 #ifndef NDEBUG
431  found_var_name = true;
432 #endif
433  const std::size_t assignment_i =
434  subexpression.find(":", varname_i+1);
435 
436  libmesh_assert_not_equal_to(assignment_i, std::string::npos);
437 
438  libmesh_assert_equal_to(subexpression[assignment_i+1], '=');
439  for (unsigned int i = varname_i+1; i != assignment_i; ++i)
440  libmesh_assert_equal_to(subexpression[i], ' ');
441 
442  std::size_t end_assignment_i =
443  subexpression.find(";", assignment_i+1);
444 
445  libmesh_assert_not_equal_to(end_assignment_i, std::string::npos);
446 
447  std::ostringstream new_subexpression;
448  new_subexpression << subexpression.substr(0, assignment_i+2)
449  << std::setprecision(std::numeric_limits<Output>::digits10+2)
450 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
451  << '(' << newval.real() << '+'
452  << newval.imag() << 'i' << ')'
453 #else
454  << newval
455 #endif
456  << subexpression.substr(end_assignment_i,
457  std::string::npos);
458  _subexpressions[s] = new_subexpression.str();
459  }
460 
461  libmesh_assert(found_var_name);
462 
463  std::string new_expression;
464 
465  for (std::size_t s=0; s != _subexpressions.size(); ++s)
466  {
467  new_expression += '{';
468  new_expression += _subexpressions[s];
469  new_expression += '}';
470  }
471 
472  this->partial_reparse(new_expression);
473 }
474 
475 
476 template <typename Output, typename OutputGradient>
477 inline
478 void
480 {
482  _subexpressions.clear();
483  parsers.clear();
484 
485  size_t nextstart = 0, end = 0;
486 
487  while (end != std::string::npos)
488  {
489  // If we're past the end of the string, we can't make any more
490  // subparsers
491  if (nextstart >= expression.size())
492  break;
493 
494  // If we're at the start of a brace delimited section, then we
495  // parse just that section:
496  if (expression[nextstart] == '{')
497  {
498  nextstart++;
499  end = expression.find('}', nextstart);
500  }
501  // otherwise we parse the whole thing
502  else
503  end = std::string::npos;
504 
505  // We either want the whole end of the string (end == npos) or
506  // a substring in the middle.
507  _subexpressions.push_back
508  (expression.substr(nextstart, (end == std::string::npos) ?
509  std::string::npos : end - nextstart));
510 
511  // fparser can crash on empty expressions
512  if (_subexpressions.back().empty())
513  libmesh_error_msg("ERROR: FunctionParser is unable to parse empty expression.\n");
514 
515  // Parse (and optimize if possible) the subexpression.
516  // Add some basic constants, to Real precision.
517  FunctionParserADBase<Output> fp;
518  fp.AddConstant("NaN", std::numeric_limits<Real>::quiet_NaN());
519  fp.AddConstant("pi", std::acos(Real(-1)));
520  fp.AddConstant("e", std::exp(Real(1)));
521  if (fp.Parse(_subexpressions.back(), variables) != -1) // -1 for success
522  libmesh_error_msg
523  ("ERROR: FunctionParser is unable to parse expression: "
524  << _subexpressions.back() << '\n' << fp.ErrorMsg());
525 
526  // use of derivatives is optional. suppress error output on the console
527  // use the has_derivatives() method to check if AutoDiff was successful.
528  // also enable immediate optimization
529  fp.SetADFlags(FunctionParserADBase<Output>::ADSilenceErrors |
530  FunctionParserADBase<Output>::ADAutoOptimize);
531 
532  // optimize original function
533  fp.Optimize();
534  parsers.push_back(fp);
535 
536  // generate derivatives through automatic differentiation
537  FunctionParserADBase<Output> dx_fp(fp);
538  if (dx_fp.AutoDiff("x") != -1) // -1 for success
539  _valid_derivatives = false;
540  dx_parsers.push_back(dx_fp);
541 #if LIBMESH_DIM > 1
542  FunctionParserADBase<Output> dy_fp(fp);
543  if (dy_fp.AutoDiff("y") != -1) // -1 for success
544  _valid_derivatives = false;
545  dy_parsers.push_back(dy_fp);
546 #endif
547 #if LIBMESH_DIM > 2
548  FunctionParserADBase<Output> dz_fp(fp);
549  if (dz_fp.AutoDiff("z") != -1) // -1 for success
550  _valid_derivatives = false;
551  dz_parsers.push_back(dz_fp);
552 #endif
553  FunctionParserADBase<Output> dt_fp(fp);
554  if (dt_fp.AutoDiff("t") != -1) // -1 for success
555  _valid_derivatives = false;
556  dt_parsers.push_back(dt_fp);
557 
558  // If at end, use nextstart=maxSize. Else start at next
559  // character.
560  nextstart = (end == std::string::npos) ?
561  std::string::npos : end + 1;
562  }
563 }
564 
565 
566 template <typename Output, typename OutputGradient>
567 inline
568 std::size_t
570  const std::string & expr) const
571 {
572  const std::size_t namesize = varname.size();
573  std::size_t varname_i = expr.find(varname);
574 
575  while ((varname_i != std::string::npos) &&
576  (((varname_i > 0) &&
577  (std::isalnum(expr[varname_i-1]) ||
578  (expr[varname_i-1] == '_'))) ||
579  ((varname_i+namesize < expr.size()) &&
580  (std::isalnum(expr[varname_i+namesize]) ||
581  (expr[varname_i+namesize] == '_')))))
582  {
583  varname_i = expr.find(varname, varname_i+1);
584  }
585 
586  return varname_i;
587 }
588 template <typename Output, typename OutputGradient>
589 inline
590 bool
592 {
593  bool is_time_dependent = false;
594 
595  // By definition, time is "t" for FunctionBase-based objects, so we just need to
596  // see if this expression has the variable "t" in it.
597  if (this->find_name(std::string("t"), expression) != std::string::npos)
598  is_time_dependent = true;
599 
600  return is_time_dependent;
601 }
602 
603 // Set the _spacetime argument vector
604 template <typename Output, typename OutputGradient>
605 inline
606 void
608  const Real time)
609 {
610  _spacetime[0] = p(0);
611 #if LIBMESH_DIM > 1
612  _spacetime[1] = p(1);
613 #endif
614 #if LIBMESH_DIM > 2
615  _spacetime[2] = p(2);
616 #endif
617  _spacetime[LIBMESH_DIM] = time;
618 
619  // The remaining locations in _spacetime are currently fixed at construction
620  // but could potentially be made dynamic
621 }
622 
623 // Evaluate the ith FunctionParser and check the result
624 template <typename Output, typename OutputGradient>
625 inline
626 Output
627 ParsedFunction<Output,OutputGradient>::eval (FunctionParserADBase<Output> & parser,
628  const std::string & libmesh_dbg_var(function_name),
629  unsigned int libmesh_dbg_var(component_idx)) const
630 {
631 #ifndef NDEBUG
632  Output result = parser.Eval(&_spacetime[0]);
633  int error_code = parser.EvalError();
634  if (error_code)
635  {
636  libMesh::err << "ERROR: FunctionParser is unable to evaluate component "
637  << component_idx
638  << " of expression '"
639  << function_name
640  << "' with arguments:\n";
641  for (std::size_t j=0; j<_spacetime.size(); ++j)
642  libMesh::err << '\t' << _spacetime[j] << '\n';
643  libMesh::err << '\n';
644 
645  // Currently no API to report error messages, we'll do it manually
646  std::string error_message = "Reason: ";
647 
648  switch (error_code)
649  {
650  case 1:
651  error_message += "Division by zero";
652  break;
653  case 2:
654  error_message += "Square Root error (negative value)";
655  break;
656  case 3:
657  error_message += "Log error (negative value)";
658  break;
659  case 4:
660  error_message += "Trigonometric error (asin or acos of illegal value)";
661  break;
662  case 5:
663  error_message += "Maximum recursion level reached";
664  break;
665  default:
666  error_message += "Unknown";
667  break;
668  }
669  libmesh_error_msg(error_message);
670  }
671 
672  return result;
673 #else
674  return parser.Eval(&_spacetime[0]);
675 #endif
676 }
677 
678 } // namespace libMesh
679 
680 
681 #else // !LIBMESH_HAVE_FPARSER
682 
683 
684 namespace libMesh {
685 
686 
687 template <typename Output=Number>
688 class ParsedFunction : public FunctionBase<Output>
689 {
690 public:
691  ParsedFunction (const std::string & /* expression */,
692  const std::vector<std::string> * = libmesh_nullptr,
693  const std::vector<Output> * = libmesh_nullptr) : _dummy(0)
694  {
695  libmesh_not_implemented();
696  }
697 
698  virtual Output operator() (const Point &,
699  const Real /* time */ = 0)
700  { return 0.; }
701 
702  virtual void operator() (const Point &,
703  const Real /* time */,
704  DenseVector<Output> & /* output */) {}
705 
706  virtual void init() {}
707  virtual void clear() {}
708  virtual Output & getVarAddress(const std::string & /*variable_name*/) { return _dummy; }
710  {
712  }
713 private:
714  Output _dummy;
715 };
716 
717 
718 
719 } // namespace libMesh
720 
721 
722 #endif // LIBMESH_HAVE_FPARSER
723 
724 #endif // LIBMESH_PARSED_FUNCTION_H
Output get_inline_value(const std::string &inline_var_name) const
virtual Output dot(const Point &p, const Real time=0)
OStreamProxy err
std::vector< FunctionParserADBase< Output > > dt_parsers
virtual Output & getVarAddress(const std::string &variable_name)
ParsedFunction(const std::string &expression, const std::vector< std::string > *additional_vars=libmesh_nullptr, const std::vector< Output > *initial_vals=libmesh_nullptr)
virtual void clear()
Clears the function.
A Function generated (via FParser) by parsing a mathematical expression.
std::vector< Output > _spacetime
Output eval(FunctionParserADBase< Output > &parser, const std::string &libmesh_dbg_var(function_name), unsigned int libmesh_dbg_var(component_idx)) const
Evaluate the ith FunctionParser and check the result.
std::vector< std::string > _additional_vars
std::vector< FunctionParserADBase< Output > > parsers
void set_spacetime(const Point &p, const Real time=0)
Set the _spacetime argument vector.
bool _initialized
When init() was called so that everything is ready for calls to operator() (...), then this bool is t...
const class libmesh_nullptr_t libmesh_nullptr
std::vector< std::string > _subexpressions
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
virtual Output & getVarAddress(const std::string &)
The libMesh namespace provides an interface to certain functionality in the library.
bool is_time_dependent() const
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
std::vector< FunctionParserADBase< Output > > dz_parsers
std::vector< FunctionParserADBase< Output > > dx_parsers
std::vector< FunctionParserADBase< Output > > dy_parsers
ParsedFunction(const std::string &, const std::vector< std::string > *=libmesh_nullptr, const std::vector< Output > *=libmesh_nullptr)
void reparse(const std::string &expression)
Re-parse with new expression.
void set_inline_value(const std::string &inline_var_name, Output newval)
Changes the value of an inline variable.
PetscErrorCode Vec Mat libmesh_dbg_var(j)
bool _is_time_dependent
Cache whether or not this function is actually time-dependent.
bool expression_is_time_dependent(const std::string &expression) const
virtual Output component(unsigned int i, const Point &p, Real time) libmesh_override
virtual bool has_derivatives()
Query if the automatic derivative generation was successful.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void partial_reparse(const std::string &expression)
Re-parse with minor changes to expression.
std::size_t find_name(const std::string &varname, const std::string &expr) const
Helper function for parsing out variable names.
virtual UniquePtr< FunctionBase< Output > > clone() const libmesh_override
virtual Output operator()(const Point &p, const Real time=0) libmesh_override
Defines a dense vector for use in Finite Element-type computations.
virtual UniquePtr< FunctionBase< Output > > clone() const
This is the base class for functor-like classes.
std::vector< Output > _initial_vals
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
const std::string & expression()
virtual void init()
The actual initialization process.
virtual OutputGradient gradient(const Point &p, const Real time=0)