20 #ifndef LIBMESH_PARSED_FEM_FUNCTION_H 21 #define LIBMESH_PARSED_FEM_FUNCTION_H 23 #include "libmesh/libmesh_config.h" 26 #include "libmesh/fem_function_base.h" 27 #include "libmesh/int_range.h" 28 #include "libmesh/point.h" 29 #include "libmesh/system.h" 31 #ifdef LIBMESH_HAVE_FPARSER 33 #include "libmesh/fparser.hh" 56 template <
typename Output=Number>
67 const std::vector<std::string> * additional_vars=
nullptr,
68 const std::vector<Output> * initial_vals=
nullptr);
91 virtual std::unique_ptr<FEMFunctionBase<Output>>
clone ()
const override;
95 const Real time = 0.)
override;
105 Real time=0.)
override;
137 std::size_t
find_name (std::string_view varname,
138 std::string_view expr)
const;
146 #ifdef LIBMESH_HAVE_FPARSER 147 inline Output
eval(FunctionParserBase<Output> & parser,
148 std::string_view libmesh_dbg_var(function_name),
149 unsigned int libmesh_dbg_var(component_idx))
const;
150 #else // LIBMESH_HAVE_FPARSER 151 inline Output
eval(
char & libmesh_dbg_var(parser),
152 std::string_view libmesh_dbg_var(function_name),
153 unsigned int libmesh_dbg_var(component_idx))
const;
165 #ifdef LIBMESH_HAVE_FPARSER 166 std::vector<std::unique_ptr<FunctionParserBase<Output>>>
parsers;
180 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 184 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 195 template <
typename Output>
198 std::string expression,
199 const std::vector<std::string> * additional_vars,
200 const std::vector<Output> * initial_vals) :
204 _n_requested_vars(0),
205 _n_requested_grad_components(0),
206 _n_requested_hess_components(0),
207 _requested_normals(false),
208 _need_var(_n_vars, false),
209 _need_var_grad(_n_vars*LIBMESH_DIM, false),
210 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
211 _need_var_hess(_n_vars*LIBMESH_DIM*LIBMESH_DIM, false),
213 _additional_vars (additional_vars ? *additional_vars :
std::vector<
std::string>()),
214 _initial_vals (initial_vals ? *initial_vals :
std::vector<Output>())
216 this->reparse(std::move(expression));
220 template <
typename Output>
225 _expression(other._expression),
226 _subexpressions(other._subexpressions),
227 _n_vars(other._n_vars),
228 _n_requested_vars(other._n_requested_vars),
229 _n_requested_grad_components(other._n_requested_grad_components),
230 _n_requested_hess_components(other._n_requested_hess_components),
231 _requested_normals(other._requested_normals),
232 _spacetime(other._spacetime),
233 _need_var(other._need_var),
234 _need_var_grad(other._need_var_grad),
235 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
236 _need_var_hess(other._need_var_hess),
238 variables(other.variables),
239 _additional_vars(other._additional_vars),
240 _initial_vals(other._initial_vals)
242 this->reparse(_expression);
246 template <
typename Output>
256 std::swap(tmp, *
this);
261 template <
typename Output>
275 for (
unsigned int v=0; v != _n_vars; ++v)
277 const std::string & varname = _sys.variable_name(v);
278 std::size_t varname_i = find_name(varname, expression);
282 if (varname_i == std::string::npos)
287 variables += varname;
291 for (
unsigned int v=0; v != _n_vars; ++v)
293 const std::string & varname = _sys.variable_name(v);
295 for (
unsigned int d=0; d != LIBMESH_DIM; ++d)
297 std::string gradname = std::string(
"grad_") +
298 "xyz"[d] +
'_' + varname;
299 std::size_t gradname_i = find_name(gradname, expression);
303 if (gradname_i == std::string::npos)
306 _need_var_grad[v*LIBMESH_DIM+d] =
true;
308 variables += gradname;
309 _n_requested_grad_components++;
313 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 314 for (
unsigned int v=0; v != _n_vars; ++v)
316 const std::string & varname = _sys.variable_name(v);
318 for (
unsigned int d1=0; d1 != LIBMESH_DIM; ++d1)
319 for (
unsigned int d2=0; d2 != LIBMESH_DIM; ++d2)
321 std::string hessname = std::string(
"hess_") +
322 "xyz"[d1] +
"xyz"[d2] +
'_' + varname;
323 std::size_t hessname_i = find_name(hessname, expression);
327 if (hessname_i == std::string::npos)
330 _need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+d1*LIBMESH_DIM+d2]
333 variables += hessname;
334 _n_requested_hess_components++;
337 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 340 std::size_t nx_i = find_name(
"n_x", expression);
341 std::size_t ny_i = find_name(
"n_y", expression);
342 std::size_t nz_i = find_name(
"n_z", expression);
346 if (nx_i != std::string::npos ||
347 ny_i != std::string::npos ||
348 nz_i != std::string::npos)
350 _requested_normals =
true;
367 (LIBMESH_DIM + 1 + _n_requested_vars +
368 _n_requested_grad_components + _n_requested_hess_components +
369 (_requested_normals ? LIBMESH_DIM : 0) +
370 _additional_vars.size());
375 unsigned int offset = LIBMESH_DIM + 1 + _n_requested_vars +
376 _n_requested_grad_components + _n_requested_hess_components;
380 variables +=
"," + _additional_vars[i];
383 _spacetime[offset + i] =
384 (i < _initial_vals.size()) ? _initial_vals[i] : 0;
387 this->partial_reparse(std::move(expression));
390 template <
typename Output>
395 for (
unsigned int v=0; v != _n_vars; ++v)
399 if (_n_requested_vars)
401 if (_n_requested_grad_components)
403 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 404 if (_n_requested_hess_components)
406 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 409 if (_requested_normals)
422 template <
typename Output>
424 std::unique_ptr<FEMFunctionBase<Output>>
427 return std::make_unique<ParsedFEMFunction>
428 (_sys, _expression, &_additional_vars, &_initial_vals);
431 template <
typename Output>
438 eval_args(c, p, time);
440 return eval(*parsers[0],
"f", 0);
445 template <
typename Output>
453 eval_args(c, p, time);
455 unsigned int size = output.
size();
457 libmesh_assert_equal_to (size, parsers.size());
459 for (
unsigned int i=0; i != size; ++i)
460 output(i) = eval(*parsers[i],
"f", i);
464 template <
typename Output>
472 eval_args(c, p, time);
474 libmesh_assert_less (i, parsers.size());
475 return eval(*parsers[i],
"f", i);
478 template <
typename Output>
483 libmesh_assert_greater (_subexpressions.size(), 0);
486 bool found_var_name =
false;
488 Output old_var_value(0.);
490 for (
const std::string & subexpression : _subexpressions)
492 const std::size_t varname_i =
493 find_name(inline_var_name, subexpression);
494 if (varname_i == std::string::npos)
497 const std::size_t assignment_i =
498 subexpression.find(
":", varname_i+1);
500 libmesh_assert_not_equal_to(assignment_i, std::string::npos);
502 libmesh_assert_equal_to(subexpression[assignment_i+1],
'=');
503 for (std::size_t i = varname_i+1; i != assignment_i; ++i)
504 libmesh_assert_equal_to(subexpression[i],
' ');
506 std::size_t end_assignment_i =
507 subexpression.find(
";", assignment_i+1);
509 libmesh_assert_not_equal_to(end_assignment_i, std::string::npos);
511 std::string new_subexpression =
512 subexpression.substr(0, end_assignment_i+1).
513 append(inline_var_name);
515 #ifdef LIBMESH_HAVE_FPARSER 518 auto fp = std::make_unique<FunctionParserBase<Output>>();
519 fp->AddConstant(
"NaN", std::numeric_limits<Real>::quiet_NaN());
520 fp->AddConstant(
"pi", std::acos(
Real(-1)));
521 fp->AddConstant(
"e", std::exp(
Real(1)));
523 (fp->Parse(new_subexpression, variables) != -1,
524 "ERROR: FunctionParser is unable to parse modified expression: " 525 << new_subexpression <<
'\n' << fp->ErrorMsg());
527 Output new_var_value = this->eval(*fp, new_subexpression, 0);
529 return new_var_value;
533 libmesh_assert_equal_to(old_var_value, new_var_value);
537 old_var_value = new_var_value;
538 found_var_name =
true;
543 libmesh_error_msg(
"ERROR: This functionality requires fparser!");
548 return old_var_value;
551 template <
typename Output>
557 libmesh_assert_greater (_subexpressions.size(), 0);
560 bool found_var_name =
false;
564 const std::string & subexpression = _subexpressions[s];
565 const std::size_t varname_i =
566 find_name(inline_var_name, subexpression);
567 if (varname_i == std::string::npos)
571 found_var_name =
true;
574 const std::size_t assignment_i =
575 subexpression.find(
":", varname_i+1);
577 libmesh_assert_not_equal_to(assignment_i, std::string::npos);
579 libmesh_assert_equal_to(subexpression[assignment_i+1],
'=');
580 for (std::size_t i = varname_i+1; i != assignment_i; ++i)
581 libmesh_assert_equal_to(subexpression[i],
' ');
583 std::size_t end_assignment_i =
584 subexpression.find(
";", assignment_i+1);
586 libmesh_assert_not_equal_to(end_assignment_i, std::string::npos);
588 std::ostringstream new_subexpression;
589 new_subexpression << subexpression.substr(0, assignment_i+2)
590 << std::setprecision(std::numeric_limits<Output>::digits10+2)
591 #ifdef LIBMESH_USE_COMPLEX_NUMBERS 592 <<
'(' << newval.real() <<
'+' 593 << newval.imag() <<
'i' <<
')' 597 << subexpression.substr(end_assignment_i,
599 _subexpressions[s] = new_subexpression.str();
604 std::string new_expression;
606 for (
const auto & subexpression : _subexpressions)
608 new_expression +=
'{';
609 new_expression += subexpression;
610 new_expression +=
'}';
613 this->partial_reparse(new_expression);
617 template <
typename Output>
622 _expression = std::move(expression);
623 _subexpressions.clear();
626 size_t nextstart = 0, end = 0;
628 while (end != std::string::npos)
632 if (nextstart >= _expression.size())
637 if (_expression[nextstart] ==
'{')
640 end = _expression.find(
'}', nextstart);
644 end = std::string::npos;
648 _subexpressions.push_back
649 (_expression.substr(nextstart, (end == std::string::npos) ?
650 std::string::npos : end - nextstart));
653 libmesh_error_msg_if(_subexpressions.back().empty(),
654 "ERROR: FunctionParser is unable to parse empty expression.\n");
657 #ifdef LIBMESH_HAVE_FPARSER 660 auto fp = std::make_unique<FunctionParserBase<Output>>();
661 fp->AddConstant(
"NaN", std::numeric_limits<Real>::quiet_NaN());
662 fp->AddConstant(
"pi", std::acos(
Real(-1)));
663 fp->AddConstant(
"e", std::exp(
Real(1)));
665 (fp->Parse(_subexpressions.back(), variables) != -1,
666 "ERROR: FunctionParser is unable to parse expression: " 667 << _subexpressions.back() <<
'\n' << fp->ErrorMsg());
669 parsers.push_back(std::move(fp));
671 libmesh_error_msg(
"ERROR: This functionality requires fparser!");
676 nextstart = (end == std::string::npos) ?
677 std::string::npos : end + 1;
681 template <
typename Output>
685 std::string_view expr)
const 687 const std::size_t namesize = varname.size();
688 std::size_t varname_i = expr.find(varname);
690 while ((varname_i != std::string::npos) &&
692 (std::isalnum(expr[varname_i-1]) ||
693 (expr[varname_i-1] ==
'_'))) ||
694 ((varname_i+namesize < expr.size()) &&
695 (std::isalnum(expr[varname_i+namesize]) ||
696 (expr[varname_i+namesize] ==
'_')))))
698 varname_i = expr.find(varname, varname_i+1);
707 template <
typename Output>
714 _spacetime[0] = p(0);
716 _spacetime[1] = p(1);
719 _spacetime[2] = p(2);
721 _spacetime[LIBMESH_DIM] = time;
723 unsigned int request_index = 0;
724 for (
unsigned int v=0; v != _n_vars; ++v)
729 c.
point_value(v, p, _spacetime[LIBMESH_DIM+1+request_index]);
733 if (_n_requested_grad_components)
734 for (
unsigned int v=0; v != _n_vars; ++v)
736 if (!_need_var_grad[v*LIBMESH_DIM]
738 && !_need_var_grad[v*LIBMESH_DIM+1]
740 && !_need_var_grad[v*LIBMESH_DIM+2]
749 for (
unsigned int d=0; d != LIBMESH_DIM; ++d)
751 if (!_need_var_grad[v*LIBMESH_DIM+d])
754 _spacetime[LIBMESH_DIM+1+request_index] = g(d);
759 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 760 if (_n_requested_hess_components)
761 for (
unsigned int v=0; v != _n_vars; ++v)
763 if (!_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM]
765 && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+1]
766 && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+2]
767 && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+3]
769 && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+4]
770 && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+5]
771 && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+6]
772 && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+7]
773 && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+8]
782 for (
unsigned int d1=0; d1 != LIBMESH_DIM; ++d1)
783 for (
unsigned int d2=0; d2 != LIBMESH_DIM; ++d2)
785 if (!_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+d1*LIBMESH_DIM+d2])
788 _spacetime[LIBMESH_DIM+1+request_index] = h(d1,d2);
792 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 794 if (_requested_normals)
799 const std::vector<Point> & normals = side_fe->
get_normals();
801 const std::vector<Point> & xyz = side_fe->
get_xyz();
803 libmesh_assert_equal_to(normals.size(), xyz.size());
807 bool at_quadrature_point =
false;
813 const Point & n = normals[qp];
814 for (
unsigned int d=0; d != LIBMESH_DIM; ++d)
816 _spacetime[LIBMESH_DIM+1+request_index] = n(d);
820 at_quadrature_point =
true;
835 #ifdef LIBMESH_HAVE_FPARSER 836 template <
typename Output>
840 std::string_view libmesh_dbg_var(function_name),
841 unsigned int libmesh_dbg_var(component_idx))
const 844 Output result = parser.Eval(_spacetime.data());
845 int error_code = parser.EvalError();
848 libMesh::err <<
"ERROR: FunctionParser is unable to evaluate component " 854 if (parsers[i].
get() == &parser)
856 << _subexpressions[i];
859 for (
const auto & st : _spacetime)
864 std::string error_message =
"Reason: ";
869 error_message +=
"Division by zero";
872 error_message +=
"Square Root error (negative value)";
875 error_message +=
"Log error (negative value)";
878 error_message +=
"Trigonometric error (asin or acos of illegal value)";
881 error_message +=
"Maximum recursion level reached";
884 error_message +=
"Unknown";
887 libmesh_error_msg(error_message);
892 return parser.Eval(_spacetime.data());
895 #else // LIBMESH_HAVE_FPARSER 896 template <
typename Output>
903 libmesh_error_msg(
"ERROR: This functionality requires fparser!");
911 #endif // LIBMESH_PARSED_FEM_FUNCTION_H
unsigned int _n_requested_vars
std::vector< bool > _need_var_grad
virtual void init_context(const FEMContext &c) override
Prepares a context object for use.
std::vector< Output > _initial_vals
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
std::vector< std::string > _additional_vars
std::size_t find_name(std::string_view varname, std::string_view expr) const
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
std::vector< std::string > _subexpressions
void partial_reparse(std::string expression)
std::vector< Output > _spacetime
The libMesh namespace provides an interface to certain functionality in the library.
virtual std::unique_ptr< FEMFunctionBase< Output > > clone() const override
Number point_value(unsigned int var, const Point &p) const
unsigned int _n_requested_grad_components
virtual Output operator()(const FEMContext &c, const Point &p, const Real time=0.) override
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Output get_inline_value(std::string_view inline_var_name) const
ParsedFEMFunction provides support for FParser-based parsed functions in FEMSystem.
Manages consistently variables, degrees of freedom, and coefficient vectors.
void reparse(std::string expression)
Re-parse with new expression.
virtual ~ParsedFEMFunction()=default
std::vector< std::unique_ptr< FunctionParserBase< Output > > > parsers
Tensor point_hessian(unsigned int var, const Point &p) const
This class provides all data required for a physics package (e.g.
virtual_for_inffe const std::vector< Point > & get_xyz() const
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
ParsedFEMFunction(const System &sys, std::string expression, const std::vector< std::string > *additional_vars=nullptr, const std::vector< Output > *initial_vals=nullptr)
Constructor.
ParsedFEMFunction & operator=(const ParsedFEMFunction &)
Constructors.
std::vector< bool > _need_var_hess
void set_inline_value(std::string_view inline_var_name, Output newval)
Changes the value of an inline variable.
void eval_args(const FEMContext &c, const Point &p, const Real time)
unsigned int _n_requested_hess_components
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< bool > _need_var
virtual_for_inffe const std::vector< Point > & get_normals() const
virtual unsigned int size() const override final
virtual Output component(const FEMContext &c, unsigned int i, const Point &p, Real time=0.) override
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
std::vector< char * > parsers
Gradient point_gradient(unsigned int var, const Point &p) const
FEMFunctionBase is a base class from which users can derive in order to define "function-like" object...
A Point defines a location in LIBMESH_DIM dimensional Real space.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
This class forms the foundation from which generic finite elements may be derived.
const std::string & expression()
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
Output eval(FunctionParserBase< Output > &parser, std::string_view libmesh_dbg_var(function_name), unsigned int libmesh_dbg_var(component_idx)) const
const std::vector< std::vector< OutputShape > > & get_phi() const