parsed_function.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 #include "libmesh/auto_ptr.h" // libmesh_make_unique
31 
32 // FParser includes
33 #include "libmesh/fparser_ad.hh"
34 
35 // C++ includes
36 #include <algorithm> // std::find
37 #include <cmath>
38 #include <cmath>
39 #include <cstddef>
40 #include <iomanip>
41 #include <limits>
42 #include <sstream>
43 #include <string>
44 #include <vector>
45 
46 namespace libMesh
47 {
48 
58 template <typename Output=Number, typename OutputGradient=Gradient>
59 class ParsedFunction : public FunctionBase<Output>
60 {
61 public:
62  explicit
63  ParsedFunction (const std::string & expression,
64  const std::vector<std::string> * additional_vars=nullptr,
65  const std::vector<Output> * initial_vals=nullptr);
66 
72  ParsedFunction & operator= (const ParsedFunction &) = delete;
73 
88  ParsedFunction (const ParsedFunction &) = default;
89  ParsedFunction (ParsedFunction &&) = default;
91  virtual ~ParsedFunction () = default;
92 
96  void reparse (const std::string & expression);
97 
98  virtual Output operator() (const Point & p,
99  const Real time = 0) override;
100 
104  virtual bool has_derivatives() { return _valid_derivatives; }
105 
106  virtual Output dot(const Point & p,
107  const Real time = 0);
108 
109  virtual OutputGradient gradient(const Point & p,
110  const Real time = 0);
111 
112  virtual void operator() (const Point & p,
113  const Real time,
114  DenseVector<Output> & output) override;
115 
116  virtual Output component (unsigned int i,
117  const Point & p,
118  Real time) override;
119 
120  const std::string & expression() { return _expression; }
121 
125  virtual Output & getVarAddress(const std::string & variable_name);
126 
127  virtual std::unique_ptr<FunctionBase<Output>> clone() const override;
128 
137  Output get_inline_value(const std::string & inline_var_name) const;
138 
149  void set_inline_value(const std::string & inline_var_name,
150  Output newval);
151 
152 protected:
156  void partial_reparse (const std::string & expression);
157 
161  std::size_t find_name (const std::string & varname,
162  const std::string & expr) const;
163 
167  bool expression_is_time_dependent( const std::string & expression ) const;
168 
169 private:
173  void set_spacetime(const Point & p,
174  const Real time = 0);
175 
179  inline Output eval(FunctionParserADBase<Output> & parser,
180  const std::string & libmesh_dbg_var(function_name),
181  unsigned int libmesh_dbg_var(component_idx)) const;
182 
183  std::string _expression;
184  std::vector<std::string> _subexpressions;
185  std::vector<FunctionParserADBase<Output>> parsers;
186  std::vector<Output> _spacetime;
187 
188  // derivative functions
189  std::vector<FunctionParserADBase<Output>> dx_parsers;
190 #if LIBMESH_DIM > 1
191  std::vector<FunctionParserADBase<Output>> dy_parsers;
192 #endif
193 #if LIBMESH_DIM > 2
194  std::vector<FunctionParserADBase<Output>> dz_parsers;
195 #endif
196  std::vector<FunctionParserADBase<Output>> dt_parsers;
198 
199  // Variables/values that can be parsed and handled by the function parser
200  std::string variables;
201  std::vector<std::string> _additional_vars;
202  std::vector<Output> _initial_vals;
203 };
204 
205 
206 /*----------------------- Inline functions ----------------------------------*/
207 
208 template <typename Output, typename OutputGradient>
209 inline
211  const std::vector<std::string> * additional_vars,
212  const std::vector<Output> * initial_vals) :
213  _expression (), // overridden by parse()
214  // Size the spacetime vector to account for space, time, and any additional
215  // variables passed
216  _spacetime (LIBMESH_DIM+1 + (additional_vars ? additional_vars->size() : 0)),
217  _valid_derivatives (true),
218  _additional_vars (additional_vars ? *additional_vars : std::vector<std::string>()),
219  _initial_vals (initial_vals ? *initial_vals : std::vector<Output>())
220 {
221  // time-dependence established in reparse function
222  this->reparse(expression);
223  this->_initialized = true;
224 }
225 
226 
227 template <typename Output, typename OutputGradient>
228 inline
229 void
230 ParsedFunction<Output,OutputGradient>::reparse (const std::string & expression)
231 {
232  variables = "x";
233 #if LIBMESH_DIM > 1
234  variables += ",y";
235 #endif
236 #if LIBMESH_DIM > 2
237  variables += ",z";
238 #endif
239  variables += ",t";
240 
241  // If additional vars were passed, append them to the string
242  // that we send to the function parser. Also add them to the
243  // end of our spacetime vector
244  for (std::size_t i=0; i < _additional_vars.size(); ++i)
245  {
246  variables += "," + _additional_vars[i];
247  // Initialize extra variables to the vector passed in or zero
248  // Note: The initial_vals vector can be shorter than the additional_vars vector
249  _spacetime[LIBMESH_DIM+1 + i] =
250  (i < _initial_vals.size()) ? _initial_vals[i] : 0;
251  }
252 
253  this->partial_reparse(expression);
254 
255  this->_is_time_dependent = this->expression_is_time_dependent(expression);
256 }
257 
258 template <typename Output, typename OutputGradient>
259 inline
260 Output
262 {
263  set_spacetime(p, time);
264  return eval(parsers[0], "f", 0);
265 }
266 
267 template <typename Output, typename OutputGradient>
268 inline
269 Output
271 {
272  set_spacetime(p, time);
273  return eval(dt_parsers[0], "df/dt", 0);
274 }
275 
276 template <typename Output, typename OutputGradient>
277 inline
278 OutputGradient
280 {
281  OutputGradient grad;
282  set_spacetime(p, time);
283 
284  grad(0) = eval(dx_parsers[0], "df/dx", 0);
285 #if LIBMESH_DIM > 1
286  grad(1) = eval(dy_parsers[0], "df/dy", 0);
287 #endif
288 #if LIBMESH_DIM > 2
289  grad(2) = eval(dz_parsers[0], "df/dz", 0);
290 #endif
291 
292  return grad;
293 }
294 
295 template <typename Output, typename OutputGradient>
296 inline
297 void
299  (const Point & p,
300  const Real time,
301  DenseVector<Output> & output)
302 {
303  set_spacetime(p, time);
304 
305  unsigned int size = output.size();
306 
307  libmesh_assert_equal_to (size, parsers.size());
308 
309  // The remaining locations in _spacetime are currently fixed at construction
310  // but could potentially be made dynamic
311  for (unsigned int i=0; i != size; ++i)
312  output(i) = eval(parsers[i], "f", i);
313 }
314 
319 template <typename Output, typename OutputGradient>
320 inline
321 Output
323  const Point & p,
324  Real time)
325 {
326  set_spacetime(p, time);
327  libmesh_assert_less (i, parsers.size());
328 
329  // The remaining locations in _spacetime are currently fixed at construction
330  // but could potentially be made dynamic
331  libmesh_assert_less(i, parsers.size());
332  return eval(parsers[i], "f", i);
333 }
334 
338 template <typename Output, typename OutputGradient>
339 inline
340 Output &
341 ParsedFunction<Output,OutputGradient>::getVarAddress (const std::string & variable_name)
342 {
343  const std::vector<std::string>::iterator it =
344  std::find(_additional_vars.begin(), _additional_vars.end(), variable_name);
345 
346  if (it == _additional_vars.end())
347  libmesh_error_msg("ERROR: Requested variable not found in parsed function");
348 
349  // Iterator Arithmetic (How far from the end of the array is our target address?)
350  return _spacetime[_spacetime.size() - (_additional_vars.end() - it)];
351 }
352 
353 
354 template <typename Output, typename OutputGradient>
355 inline
356 std::unique_ptr<FunctionBase<Output>>
358 {
359  return libmesh_make_unique<ParsedFunction>(_expression,
360  &_additional_vars,
361  &_initial_vals);
362 }
363 
364 template <typename Output, typename OutputGradient>
365 inline
366 Output
367 ParsedFunction<Output,OutputGradient>::get_inline_value (const std::string & inline_var_name) const
368 {
369  libmesh_assert_greater (_subexpressions.size(), 0);
370 
371 #ifndef NDEBUG
372  bool found_var_name = false;
373 #endif
374  Output old_var_value(0.);
375 
376  for (const auto & subexpression : _subexpressions)
377  {
378  const std::size_t varname_i =
379  find_name(inline_var_name, subexpression);
380  if (varname_i == std::string::npos)
381  continue;
382 
383  const std::size_t assignment_i =
384  subexpression.find(":", varname_i+1);
385 
386  libmesh_assert_not_equal_to(assignment_i, std::string::npos);
387 
388  libmesh_assert_equal_to(subexpression[assignment_i+1], '=');
389  for (std::size_t i = varname_i+1; i != assignment_i; ++i)
390  libmesh_assert_equal_to(subexpression[i], ' ');
391 
392  std::size_t end_assignment_i =
393  subexpression.find(";", assignment_i+1);
394 
395  libmesh_assert_not_equal_to(end_assignment_i, std::string::npos);
396 
397  std::string new_subexpression =
398  subexpression.substr(0, end_assignment_i+1) +
399  inline_var_name;
400 
401 #ifdef LIBMESH_HAVE_FPARSER
402  // Parse and evaluate the new subexpression.
403  // Add the same constants as we used originally.
404  FunctionParserADBase<Output> fp;
405  fp.AddConstant("NaN", std::numeric_limits<Real>::quiet_NaN());
406  fp.AddConstant("pi", std::acos(Real(-1)));
407  fp.AddConstant("e", std::exp(Real(1)));
408  if (fp.Parse(new_subexpression, variables) != -1) // -1 for success
409  libmesh_error_msg
410  ("ERROR: FunctionParser is unable to parse modified expression: "
411  << new_subexpression << '\n' << fp.ErrorMsg());
412 
413  Output new_var_value = this->eval(fp, new_subexpression, 0);
414 #ifdef NDEBUG
415  return new_var_value;
416 #else
417  if (found_var_name)
418  {
419  libmesh_assert_equal_to(old_var_value, new_var_value);
420  }
421  else
422  {
423  old_var_value = new_var_value;
424  found_var_name = true;
425  }
426 #endif
427 
428 #else
429  libmesh_error_msg("ERROR: This functionality requires fparser!");
430 #endif
431  }
432 
433  libmesh_assert(found_var_name);
434  return old_var_value;
435 }
436 
437 
438 template <typename Output, typename OutputGradient>
439 inline
440 void
441 ParsedFunction<Output,OutputGradient>::set_inline_value (const std::string & inline_var_name,
442  Output newval)
443 {
444  libmesh_assert_greater (_subexpressions.size(), 0);
445 
446 #ifndef NDEBUG
447  bool found_var_name = false;
448 #endif
449  for (auto & subexpression : _subexpressions)
450  {
451  const std::size_t varname_i =
452  find_name(inline_var_name, subexpression);
453  if (varname_i == std::string::npos)
454  continue;
455 
456 #ifndef NDEBUG
457  found_var_name = true;
458 #endif
459  const std::size_t assignment_i =
460  subexpression.find(":", varname_i+1);
461 
462  libmesh_assert_not_equal_to(assignment_i, std::string::npos);
463 
464  libmesh_assert_equal_to(subexpression[assignment_i+1], '=');
465  for (std::size_t i = varname_i+1; i != assignment_i; ++i)
466  libmesh_assert_equal_to(subexpression[i], ' ');
467 
468  std::size_t end_assignment_i =
469  subexpression.find(";", assignment_i+1);
470 
471  libmesh_assert_not_equal_to(end_assignment_i, std::string::npos);
472 
473  std::ostringstream new_subexpression;
474  new_subexpression << subexpression.substr(0, assignment_i+2)
475  << std::setprecision(std::numeric_limits<Output>::digits10+2)
476 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
477  << '(' << newval.real() << '+'
478  << newval.imag() << 'i' << ')'
479 #else
480  << newval
481 #endif
482  << subexpression.substr(end_assignment_i,
483  std::string::npos);
484  subexpression = new_subexpression.str();
485  }
486 
487  libmesh_assert(found_var_name);
488 
489  std::string new_expression;
490 
491  for (const auto & subexpression : _subexpressions)
492  {
493  new_expression += '{';
494  new_expression += subexpression;
495  new_expression += '}';
496  }
497 
498  this->partial_reparse(new_expression);
499 }
500 
501 
502 template <typename Output, typename OutputGradient>
503 inline
504 void
506 {
507  _expression = expression;
508  _subexpressions.clear();
509  parsers.clear();
510 
511  size_t nextstart = 0, end = 0;
512 
513  while (end != std::string::npos)
514  {
515  // If we're past the end of the string, we can't make any more
516  // subparsers
517  if (nextstart >= expression.size())
518  break;
519 
520  // If we're at the start of a brace delimited section, then we
521  // parse just that section:
522  if (expression[nextstart] == '{')
523  {
524  nextstart++;
525  end = expression.find('}', nextstart);
526  }
527  // otherwise we parse the whole thing
528  else
529  end = std::string::npos;
530 
531  // We either want the whole end of the string (end == npos) or
532  // a substring in the middle.
533  _subexpressions.push_back
534  (expression.substr(nextstart, (end == std::string::npos) ?
535  std::string::npos : end - nextstart));
536 
537  // fparser can crash on empty expressions
538  if (_subexpressions.back().empty())
539  libmesh_error_msg("ERROR: FunctionParser is unable to parse empty expression.\n");
540 
541  // Parse (and optimize if possible) the subexpression.
542  // Add some basic constants, to Real precision.
543  FunctionParserADBase<Output> fp;
544  fp.AddConstant("NaN", std::numeric_limits<Real>::quiet_NaN());
545  fp.AddConstant("pi", std::acos(Real(-1)));
546  fp.AddConstant("e", std::exp(Real(1)));
547  if (fp.Parse(_subexpressions.back(), variables) != -1) // -1 for success
548  libmesh_error_msg
549  ("ERROR: FunctionParser is unable to parse expression: "
550  << _subexpressions.back() << '\n' << fp.ErrorMsg());
551 
552  // use of derivatives is optional. suppress error output on the console
553  // use the has_derivatives() method to check if AutoDiff was successful.
554  // also enable immediate optimization
555  fp.SetADFlags(FunctionParserADBase<Output>::ADSilenceErrors |
556  FunctionParserADBase<Output>::ADAutoOptimize);
557 
558  // optimize original function
559  fp.Optimize();
560  parsers.push_back(fp);
561 
562  // generate derivatives through automatic differentiation
563  FunctionParserADBase<Output> dx_fp(fp);
564  if (dx_fp.AutoDiff("x") != -1) // -1 for success
565  _valid_derivatives = false;
566  dx_parsers.push_back(dx_fp);
567 #if LIBMESH_DIM > 1
568  FunctionParserADBase<Output> dy_fp(fp);
569  if (dy_fp.AutoDiff("y") != -1) // -1 for success
570  _valid_derivatives = false;
571  dy_parsers.push_back(dy_fp);
572 #endif
573 #if LIBMESH_DIM > 2
574  FunctionParserADBase<Output> dz_fp(fp);
575  if (dz_fp.AutoDiff("z") != -1) // -1 for success
576  _valid_derivatives = false;
577  dz_parsers.push_back(dz_fp);
578 #endif
579  FunctionParserADBase<Output> dt_fp(fp);
580  if (dt_fp.AutoDiff("t") != -1) // -1 for success
581  _valid_derivatives = false;
582  dt_parsers.push_back(dt_fp);
583 
584  // If at end, use nextstart=maxSize. Else start at next
585  // character.
586  nextstart = (end == std::string::npos) ?
587  std::string::npos : end + 1;
588  }
589 }
590 
591 
592 template <typename Output, typename OutputGradient>
593 inline
594 std::size_t
596  const std::string & expr) const
597 {
598  const std::size_t namesize = varname.size();
599  std::size_t varname_i = expr.find(varname);
600 
601  while ((varname_i != std::string::npos) &&
602  (((varname_i > 0) &&
603  (std::isalnum(expr[varname_i-1]) ||
604  (expr[varname_i-1] == '_'))) ||
605  ((varname_i+namesize < expr.size()) &&
606  (std::isalnum(expr[varname_i+namesize]) ||
607  (expr[varname_i+namesize] == '_')))))
608  {
609  varname_i = expr.find(varname, varname_i+1);
610  }
611 
612  return varname_i;
613 }
614 template <typename Output, typename OutputGradient>
615 inline
616 bool
618 {
619  bool is_time_dependent = false;
620 
621  // By definition, time is "t" for FunctionBase-based objects, so we just need to
622  // see if this expression has the variable "t" in it.
623  if (this->find_name(std::string("t"), expression) != std::string::npos)
624  is_time_dependent = true;
625 
626  return is_time_dependent;
627 }
628 
629 // Set the _spacetime argument vector
630 template <typename Output, typename OutputGradient>
631 inline
632 void
634  const Real time)
635 {
636  _spacetime[0] = p(0);
637 #if LIBMESH_DIM > 1
638  _spacetime[1] = p(1);
639 #endif
640 #if LIBMESH_DIM > 2
641  _spacetime[2] = p(2);
642 #endif
643  _spacetime[LIBMESH_DIM] = time;
644 
645  // The remaining locations in _spacetime are currently fixed at construction
646  // but could potentially be made dynamic
647 }
648 
649 // Evaluate the ith FunctionParser and check the result
650 template <typename Output, typename OutputGradient>
651 inline
652 Output
653 ParsedFunction<Output,OutputGradient>::eval (FunctionParserADBase<Output> & parser,
654  const std::string & libmesh_dbg_var(function_name),
655  unsigned int libmesh_dbg_var(component_idx)) const
656 {
657 #ifndef NDEBUG
658  Output result = parser.Eval(_spacetime.data());
659  int error_code = parser.EvalError();
660  if (error_code)
661  {
662  libMesh::err << "ERROR: FunctionParser is unable to evaluate component "
663  << component_idx
664  << " of expression '"
665  << function_name
666  << "' with arguments:\n";
667  for (const auto & item : _spacetime)
668  libMesh::err << '\t' << item << '\n';
669  libMesh::err << '\n';
670 
671  // Currently no API to report error messages, we'll do it manually
672  std::string error_message = "Reason: ";
673 
674  switch (error_code)
675  {
676  case 1:
677  error_message += "Division by zero";
678  break;
679  case 2:
680  error_message += "Square Root error (negative value)";
681  break;
682  case 3:
683  error_message += "Log error (negative value)";
684  break;
685  case 4:
686  error_message += "Trigonometric error (asin or acos of illegal value)";
687  break;
688  case 5:
689  error_message += "Maximum recursion level reached";
690  break;
691  default:
692  error_message += "Unknown";
693  break;
694  }
695  libmesh_error_msg(error_message);
696  }
697 
698  return result;
699 #else
700  return parser.Eval(_spacetime.data());
701 #endif
702 }
703 
704 } // namespace libMesh
705 
706 
707 #else // !LIBMESH_HAVE_FPARSER
708 
709 
710 namespace libMesh {
711 
712 
713 template <typename Output=Number>
714 class ParsedFunction : public FunctionBase<Output>
715 {
716 public:
717  ParsedFunction (const std::string & /* expression */,
718  const std::vector<std::string> * = nullptr,
719  const std::vector<Output> * = nullptr) : _dummy(0)
720  {
721  libmesh_not_implemented();
722  }
723 
728  ParsedFunction (ParsedFunction &&) = delete;
729  ParsedFunction (const ParsedFunction &) = delete;
730  ParsedFunction & operator= (const ParsedFunction &) = delete;
732  virtual ~ParsedFunction () = default;
733 
734  virtual Output operator() (const Point &,
735  const Real /* time */ = 0)
736  { return 0.; }
737 
738  virtual void operator() (const Point &,
739  const Real /* time */,
740  DenseVector<Output> & /* output */) {}
741 
742  virtual void init() {}
743  virtual void clear() {}
744  virtual Output & getVarAddress(const std::string & /*variable_name*/) { return _dummy; }
745  virtual std::unique_ptr<FunctionBase<Output>> clone() const
746  {
747  return libmesh_make_unique<ParsedFunction<Output>>("");
748  }
749 private:
750  Output _dummy;
751 };
752 
753 
754 
755 } // namespace libMesh
756 
757 
758 #endif // LIBMESH_HAVE_FPARSER
759 
760 #endif // LIBMESH_PARSED_FUNCTION_H
virtual Output dot(const Point &p, const Real time=0)
std::vector< FunctionParserADBase< Output > > dt_parsers
virtual std::unique_ptr< FunctionBase< Output > > clone() const override
virtual Output & getVarAddress(const std::string &variable_name)
A Function defined by a std::string.
ParsedFunction(const std::string &, const std::vector< std::string > *=nullptr, const std::vector< Output > *=nullptr)
std::vector< Output > _spacetime
std::vector< std::string > _additional_vars
std::vector< FunctionParserADBase< Output > > parsers
void set_spacetime(const Point &p, const Real time=0)
std::vector< std::string > _subexpressions
IterBase * end
virtual Output & getVarAddress(const std::string &)
std::vector< FunctionParserADBase< Output > > dz_parsers
std::vector< FunctionParserADBase< Output > > dx_parsers
std::vector< FunctionParserADBase< Output > > dy_parsers
void reparse(const std::string &expression)
void set_inline_value(const std::string &inline_var_name, Output newval)
Output eval(FunctionParserADBase< Output > &parser, const std::string &libmesh_dbg_var(function_name), unsigned int libmesh_dbg_var(component_idx)) const
virtual Output component(unsigned int i, const Point &p, Real time) override
OStreamProxy err(std::cerr)
virtual ~ParsedFunction()=default
bool expression_is_time_dependent(const std::string &expression) const
std::size_t find_name(const std::string &varname, const std::string &expr) const
virtual std::unique_ptr< FunctionBase< Output > > clone() const
virtual bool has_derivatives()
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void partial_reparse(const std::string &expression)
Base class for functors that can be evaluated at a point and (optionally) time.
std::vector< Output > _initial_vals
virtual Output operator()(const Point &p, const Real time=0) override
ParsedFunction(const std::string &expression, const std::vector< std::string > *additional_vars=nullptr, const std::vector< Output > *initial_vals=nullptr)
A geometric point in (x,y,z) space.
Definition: point.h:38
Output get_inline_value(const std::string &inline_var_name) const
const std::string & expression()
virtual OutputGradient gradient(const Point &p, const Real time=0)
ParsedFunction & operator=(const ParsedFunction &)=delete