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