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 :
19 :
20 : #ifndef LIBMESH_PARSED_FEM_FUNCTION_H
21 : #define LIBMESH_PARSED_FEM_FUNCTION_H
22 :
23 : #include "libmesh/libmesh_config.h"
24 :
25 : // Local Includes
26 : #include "libmesh/fem_function_base.h"
27 : #include "libmesh/int_range.h"
28 : #include "libmesh/point.h"
29 : #include "libmesh/system.h"
30 :
31 : #ifdef LIBMESH_HAVE_FPARSER
32 : // FParser includes
33 : #include "libmesh/fparser.hh"
34 : #endif
35 :
36 : // C++ includes
37 : #include <cctype>
38 : #include <iomanip>
39 : #include <sstream>
40 : #include <string>
41 : #include <vector>
42 : #include <memory>
43 :
44 : namespace libMesh
45 : {
46 :
47 : /**
48 : * ParsedFEMFunction provides support for FParser-based parsed
49 : * functions in FEMSystem. All overridden virtual functions are
50 : * documented in fem_function_base.h.
51 : *
52 : * \author Roy Stogner
53 : * \date 2014
54 : * \brief Support for using parsed functions in FEMSystem.
55 : */
56 : template <typename Output=Number>
57 : class ParsedFEMFunction : public FEMFunctionBase<Output>
58 : {
59 : public:
60 :
61 : /**
62 : * Constructor.
63 : */
64 : explicit
65 : ParsedFEMFunction (const System & sys,
66 : std::string expression,
67 : const std::vector<std::string> * additional_vars=nullptr,
68 : const std::vector<Output> * initial_vals=nullptr);
69 :
70 : /**
71 : * Constructors
72 : * - This class contains a const reference so it can't be default
73 : * copy or move-assigned. We manually implement the former
74 : * - This class contains unique_ptrs so it can't be default copy
75 : * constructed or assigned.
76 : * - It can still be default moved and deleted.
77 : */
78 : ParsedFEMFunction & operator= (const ParsedFEMFunction &);
79 : ParsedFEMFunction & operator= (ParsedFEMFunction &&) = delete;
80 : ParsedFEMFunction (const ParsedFEMFunction &);
81 : ParsedFEMFunction (ParsedFEMFunction &&) = default;
82 0 : virtual ~ParsedFEMFunction () = default;
83 :
84 : /**
85 : * Re-parse with new expression.
86 : */
87 : void reparse (std::string expression);
88 :
89 : virtual void init_context (const FEMContext & c) override;
90 :
91 : virtual std::unique_ptr<FEMFunctionBase<Output>> clone () const override;
92 :
93 : virtual Output operator() (const FEMContext & c,
94 : const Point & p,
95 : const Real time = 0.) override;
96 :
97 : void operator() (const FEMContext & c,
98 : const Point & p,
99 : const Real time,
100 : DenseVector<Output> & output) override;
101 :
102 : virtual Output component(const FEMContext & c,
103 : unsigned int i,
104 : const Point & p,
105 : Real time=0.) override;
106 :
107 : const std::string & expression() { return _expression; }
108 :
109 : /**
110 : * \returns The value of an inline variable.
111 : *
112 : * \note Will *only* be correct if the inline variable value is
113 : * independent of input variables, if the inline variable is not
114 : * redefined within any subexpression, and if the inline variable
115 : * takes the same value within any subexpressions where it appears.
116 : */
117 : Output get_inline_value(std::string_view inline_var_name) const;
118 :
119 : /**
120 : * Changes the value of an inline variable.
121 : *
122 : * \note Forever after, the variable value will take the given
123 : * constant, independent of input variables, in every subexpression
124 : * where it is already defined.
125 : *
126 : * \note Currently only works if the inline variable is not redefined
127 : * within any one subexpression.
128 : */
129 : void set_inline_value(std::string_view inline_var_name,
130 : Output newval);
131 :
132 : protected:
133 : // Helper function for reparsing minor changes to expression
134 : void partial_reparse (std::string expression);
135 :
136 : // Helper function for parsing out variable names
137 : std::size_t find_name (std::string_view varname,
138 : std::string_view expr) const;
139 :
140 : // Helper function for evaluating function arguments
141 : void eval_args(const FEMContext & c,
142 : const Point & p,
143 : const Real time);
144 :
145 : // Evaluate the ith FunctionParser and check the result
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;
154 : #endif
155 :
156 : private:
157 : const System & _sys;
158 : std::string _expression;
159 : std::vector<std::string> _subexpressions;
160 : unsigned int _n_vars,
161 : _n_requested_vars,
162 : _n_requested_grad_components,
163 : _n_requested_hess_components;
164 : bool _requested_normals;
165 : #ifdef LIBMESH_HAVE_FPARSER
166 : std::vector<std::unique_ptr<FunctionParserBase<Output>>> parsers;
167 : #else
168 : std::vector<char*> parsers;
169 : #endif
170 : std::vector<Output> _spacetime;
171 :
172 : // Flags for which variables need to be computed
173 :
174 : // _need_var[v] is true iff value(v) is needed
175 : std::vector<bool> _need_var;
176 :
177 : // _need_var_grad[v*LIBMESH_DIM+i] is true iff grad(v,i) is needed
178 : std::vector<bool> _need_var_grad;
179 :
180 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
181 : // _need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+i*LIBMESH_DIM+j] is true
182 : // iff grad(v,i,j) is needed
183 : std::vector<bool> _need_var_hess;
184 : #endif // LIBMESH_ENABLE_SECOND_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>
196 : inline
197 0 : ParsedFEMFunction<Output>::ParsedFEMFunction (const System & sys,
198 : std::string expression,
199 : const std::vector<std::string> * additional_vars,
200 : const std::vector<Output> * initial_vals) :
201 0 : _sys(sys),
202 0 : _expression (), // overridden by parse()
203 0 : _n_vars(sys.n_vars()),
204 0 : _n_requested_vars(0),
205 0 : _n_requested_grad_components(0),
206 0 : _n_requested_hess_components(0),
207 0 : _requested_normals(false),
208 0 : _need_var(_n_vars, false),
209 0 : _need_var_grad(_n_vars*LIBMESH_DIM, false),
210 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
211 0 : _need_var_hess(_n_vars*LIBMESH_DIM*LIBMESH_DIM, false),
212 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
213 0 : _additional_vars (additional_vars ? *additional_vars : std::vector<std::string>()),
214 0 : _initial_vals (initial_vals ? *initial_vals : std::vector<Output>())
215 : {
216 0 : this->reparse(std::move(expression));
217 0 : }
218 :
219 :
220 : template <typename Output>
221 : inline
222 : ParsedFEMFunction<Output>::ParsedFEMFunction (const ParsedFEMFunction<Output> & other) :
223 : FEMFunctionBase<Output>(other),
224 : _sys(other._sys),
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),
237 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
238 : variables(other.variables),
239 : _additional_vars(other._additional_vars),
240 : _initial_vals(other._initial_vals)
241 : {
242 : this->reparse(_expression);
243 : }
244 :
245 :
246 : template <typename Output>
247 : inline
248 : ParsedFEMFunction<Output> &
249 : ParsedFEMFunction<Output>::operator= (const ParsedFEMFunction<Output> & other)
250 : {
251 : // We can only be assigned another ParsedFEMFunction defined on the same System
252 : libmesh_assert(&_sys == &other._sys);
253 :
254 : // Use copy-and-swap idiom
255 : ParsedFEMFunction<Output> tmp(other);
256 : std::swap(tmp, *this);
257 : return *this;
258 : }
259 :
260 :
261 : template <typename Output>
262 : inline
263 : void
264 0 : ParsedFEMFunction<Output>::reparse (std::string expression)
265 : {
266 0 : variables = "x";
267 : #if LIBMESH_DIM > 1
268 0 : variables += ",y";
269 : #endif
270 : #if LIBMESH_DIM > 2
271 0 : variables += ",z";
272 : #endif
273 0 : variables += ",t";
274 :
275 0 : for (unsigned int v=0; v != _n_vars; ++v)
276 : {
277 0 : const std::string & varname = _sys.variable_name(v);
278 0 : std::size_t varname_i = find_name(varname, expression);
279 :
280 : // If we didn't find our variable name then let's go to the
281 : // next.
282 0 : if (varname_i == std::string::npos)
283 0 : continue;
284 :
285 0 : _need_var[v] = true;
286 0 : variables += ',';
287 0 : variables += varname;
288 0 : _n_requested_vars++;
289 : }
290 :
291 0 : for (unsigned int v=0; v != _n_vars; ++v)
292 : {
293 0 : const std::string & varname = _sys.variable_name(v);
294 :
295 0 : for (unsigned int d=0; d != LIBMESH_DIM; ++d)
296 : {
297 0 : std::string gradname = std::string("grad_") +
298 0 : "xyz"[d] + '_' + varname;
299 0 : std::size_t gradname_i = find_name(gradname, expression);
300 :
301 : // If we didn't find that gradient component of our
302 : // variable name then let's go to the next.
303 0 : if (gradname_i == std::string::npos)
304 0 : continue;
305 :
306 0 : _need_var_grad[v*LIBMESH_DIM+d] = true;
307 0 : variables += ',';
308 0 : variables += gradname;
309 0 : _n_requested_grad_components++;
310 : }
311 : }
312 :
313 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
314 0 : for (unsigned int v=0; v != _n_vars; ++v)
315 : {
316 0 : const std::string & varname = _sys.variable_name(v);
317 :
318 0 : for (unsigned int d1=0; d1 != LIBMESH_DIM; ++d1)
319 0 : for (unsigned int d2=0; d2 != LIBMESH_DIM; ++d2)
320 : {
321 0 : std::string hessname = std::string("hess_") +
322 0 : "xyz"[d1] + "xyz"[d2] + '_' + varname;
323 0 : std::size_t hessname_i = find_name(hessname, expression);
324 :
325 : // If we didn't find that hessian component of our
326 : // variable name then let's go to the next.
327 0 : if (hessname_i == std::string::npos)
328 0 : continue;
329 :
330 0 : _need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+d1*LIBMESH_DIM+d2]
331 0 : = true;
332 0 : variables += ',';
333 0 : variables += hessname;
334 0 : _n_requested_hess_components++;
335 : }
336 : }
337 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
338 :
339 : {
340 0 : std::size_t nx_i = find_name("n_x", expression);
341 0 : std::size_t ny_i = find_name("n_y", expression);
342 0 : std::size_t nz_i = find_name("n_z", expression);
343 :
344 : // If we found any requests for normal components then we'll
345 : // compute normals
346 0 : if (nx_i != std::string::npos ||
347 0 : ny_i != std::string::npos ||
348 : nz_i != std::string::npos)
349 : {
350 0 : _requested_normals = true;
351 0 : variables += ',';
352 0 : variables += "n_x";
353 : if (LIBMESH_DIM > 1)
354 : {
355 0 : variables += ',';
356 0 : variables += "n_y";
357 : }
358 : if (LIBMESH_DIM > 2)
359 : {
360 0 : variables += ',';
361 0 : variables += "n_z";
362 : }
363 : }
364 : }
365 :
366 0 : _spacetime.resize
367 0 : (LIBMESH_DIM + 1 + _n_requested_vars +
368 0 : _n_requested_grad_components + _n_requested_hess_components +
369 0 : (_requested_normals ? LIBMESH_DIM : 0) +
370 0 : _additional_vars.size());
371 :
372 : // If additional vars were passed, append them to the string
373 : // that we send to the function parser. Also add them to the
374 : // end of our spacetime vector
375 0 : unsigned int offset = LIBMESH_DIM + 1 + _n_requested_vars +
376 0 : _n_requested_grad_components + _n_requested_hess_components;
377 :
378 0 : for (auto i : index_range(_additional_vars))
379 : {
380 0 : variables += "," + _additional_vars[i];
381 : // Initialize extra variables to the vector passed in or zero
382 : // Note: The initial_vals vector can be shorter than the additional_vars vector
383 0 : _spacetime[offset + i] =
384 0 : (i < _initial_vals.size()) ? _initial_vals[i] : 0;
385 : }
386 :
387 0 : this->partial_reparse(std::move(expression));
388 0 : }
389 :
390 : template <typename Output>
391 : inline
392 : void
393 0 : ParsedFEMFunction<Output>::init_context (const FEMContext & c)
394 : {
395 0 : for (unsigned int v=0; v != _n_vars; ++v)
396 : {
397 : FEBase * elem_fe;
398 0 : c.get_element_fe(v, elem_fe);
399 0 : bool request_nothing = true;
400 0 : if (_n_requested_vars)
401 : {
402 0 : elem_fe->get_phi();
403 0 : request_nothing = false;
404 : }
405 0 : if (_n_requested_grad_components)
406 : {
407 0 : elem_fe->get_dphi();
408 0 : request_nothing = false;
409 : }
410 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
411 0 : if (_n_requested_hess_components)
412 : {
413 0 : elem_fe->get_d2phi();
414 0 : request_nothing = false;
415 : }
416 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
417 0 : if (request_nothing)
418 0 : elem_fe->get_nothing();
419 : }
420 :
421 0 : if (_requested_normals)
422 : {
423 : FEBase * side_fe;
424 0 : c.get_side_fe(0, side_fe);
425 :
426 0 : side_fe->get_normals();
427 :
428 : // FIXME: this is a hack to support normals at quadrature
429 : // points; we don't support normals elsewhere.
430 0 : side_fe->get_xyz();
431 : }
432 0 : }
433 :
434 : template <typename Output>
435 : inline
436 : std::unique_ptr<FEMFunctionBase<Output>>
437 0 : ParsedFEMFunction<Output>::clone () const
438 : {
439 : return std::make_unique<ParsedFEMFunction>
440 0 : (_sys, _expression, &_additional_vars, &_initial_vals);
441 : }
442 :
443 : template <typename Output>
444 : inline
445 : Output
446 0 : ParsedFEMFunction<Output>::operator() (const FEMContext & c,
447 : const Point & p,
448 : const Real time)
449 : {
450 0 : eval_args(c, p, time);
451 :
452 0 : return eval(*parsers[0], "f", 0);
453 : }
454 :
455 :
456 :
457 : template <typename Output>
458 : inline
459 : void
460 0 : ParsedFEMFunction<Output>::operator() (const FEMContext & c,
461 : const Point & p,
462 : const Real time,
463 : DenseVector<Output> & output)
464 : {
465 0 : eval_args(c, p, time);
466 :
467 0 : unsigned int size = output.size();
468 :
469 0 : libmesh_assert_equal_to (size, parsers.size());
470 :
471 0 : for (unsigned int i=0; i != size; ++i)
472 0 : output(i) = eval(*parsers[i], "f", i);
473 0 : }
474 :
475 :
476 : template <typename Output>
477 : inline
478 : Output
479 0 : ParsedFEMFunction<Output>::component (const FEMContext & c,
480 : unsigned int i,
481 : const Point & p,
482 : Real time)
483 : {
484 0 : eval_args(c, p, time);
485 :
486 0 : libmesh_assert_less (i, parsers.size());
487 0 : return eval(*parsers[i], "f", i);
488 : }
489 :
490 : template <typename Output>
491 : inline
492 : Output
493 : ParsedFEMFunction<Output>::get_inline_value(std::string_view inline_var_name) const
494 : {
495 : libmesh_assert_greater (_subexpressions.size(), 0);
496 :
497 : #ifndef NDEBUG
498 : bool found_var_name = false;
499 : #endif
500 : Output old_var_value(0.);
501 :
502 : for (const std::string & subexpression : _subexpressions)
503 : {
504 : const std::size_t varname_i =
505 : find_name(inline_var_name, subexpression);
506 : if (varname_i == std::string::npos)
507 : continue;
508 :
509 : const std::size_t assignment_i =
510 : subexpression.find(":", varname_i+1);
511 :
512 : libmesh_assert_not_equal_to(assignment_i, std::string::npos);
513 :
514 : libmesh_assert_equal_to(subexpression[assignment_i+1], '=');
515 : for (std::size_t i = varname_i+1; i != assignment_i; ++i)
516 : libmesh_assert_equal_to(subexpression[i], ' ');
517 :
518 : std::size_t end_assignment_i =
519 : subexpression.find(";", assignment_i+1);
520 :
521 : libmesh_assert_not_equal_to(end_assignment_i, std::string::npos);
522 :
523 : std::string new_subexpression =
524 : subexpression.substr(0, end_assignment_i+1).
525 : append(inline_var_name);
526 :
527 : #ifdef LIBMESH_HAVE_FPARSER
528 : // Parse and evaluate the new subexpression.
529 : // Add the same constants as we used originally.
530 : auto fp = std::make_unique<FunctionParserBase<Output>>();
531 : fp->AddConstant("NaN", std::numeric_limits<Real>::quiet_NaN());
532 : fp->AddConstant("pi", std::acos(Real(-1)));
533 : fp->AddConstant("e", std::exp(Real(1)));
534 : libmesh_error_msg_if
535 : (fp->Parse(new_subexpression, variables) != -1, // -1 for success
536 : "ERROR: FunctionParser is unable to parse modified expression: "
537 : << new_subexpression << '\n' << fp->ErrorMsg());
538 :
539 : Output new_var_value = this->eval(*fp, new_subexpression, 0);
540 : #ifdef NDEBUG
541 : return new_var_value;
542 : #else
543 : if (found_var_name)
544 : {
545 : libmesh_assert_equal_to(old_var_value, new_var_value);
546 : }
547 : else
548 : {
549 : old_var_value = new_var_value;
550 : found_var_name = true;
551 : }
552 : #endif
553 :
554 : #else
555 : libmesh_error_msg("ERROR: This functionality requires fparser!");
556 : #endif
557 : }
558 :
559 : libmesh_assert(found_var_name);
560 : return old_var_value;
561 : }
562 :
563 : template <typename Output>
564 : inline
565 : void
566 : ParsedFEMFunction<Output>::set_inline_value (std::string_view inline_var_name,
567 : Output newval)
568 : {
569 : libmesh_assert_greater (_subexpressions.size(), 0);
570 :
571 : #ifndef NDEBUG
572 : bool found_var_name = false;
573 : #endif
574 : for (auto s : index_range(_subexpressions))
575 : {
576 : const std::string & subexpression = _subexpressions[s];
577 : const std::size_t varname_i =
578 : find_name(inline_var_name, subexpression);
579 : if (varname_i == std::string::npos)
580 : continue;
581 :
582 : #ifndef NDEBUG
583 : found_var_name = true;
584 : #endif
585 :
586 : const std::size_t assignment_i =
587 : subexpression.find(":", varname_i+1);
588 :
589 : libmesh_assert_not_equal_to(assignment_i, std::string::npos);
590 :
591 : libmesh_assert_equal_to(subexpression[assignment_i+1], '=');
592 : for (std::size_t i = varname_i+1; i != assignment_i; ++i)
593 : libmesh_assert_equal_to(subexpression[i], ' ');
594 :
595 : std::size_t end_assignment_i =
596 : subexpression.find(";", assignment_i+1);
597 :
598 : libmesh_assert_not_equal_to(end_assignment_i, std::string::npos);
599 :
600 : std::ostringstream new_subexpression;
601 : new_subexpression << subexpression.substr(0, assignment_i+2)
602 : << std::setprecision(std::numeric_limits<Output>::digits10+2)
603 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
604 : << '(' << newval.real() << '+'
605 : << newval.imag() << 'i' << ')'
606 : #else
607 : << newval
608 : #endif
609 : << subexpression.substr(end_assignment_i,
610 : std::string::npos);
611 : _subexpressions[s] = new_subexpression.str();
612 : }
613 :
614 : libmesh_assert(found_var_name);
615 :
616 : std::string new_expression;
617 :
618 : for (const auto & subexpression : _subexpressions)
619 : {
620 : new_expression += '{';
621 : new_expression += subexpression;
622 : new_expression += '}';
623 : }
624 :
625 : this->partial_reparse(new_expression);
626 : }
627 :
628 :
629 : template <typename Output>
630 : inline
631 : void
632 0 : ParsedFEMFunction<Output>::partial_reparse (std::string expression)
633 : {
634 0 : _expression = std::move(expression);
635 0 : _subexpressions.clear();
636 0 : parsers.clear();
637 :
638 0 : size_t nextstart = 0, end = 0;
639 :
640 0 : while (end != std::string::npos)
641 : {
642 : // If we're past the end of the string, we can't make any more
643 : // subparsers
644 0 : if (nextstart >= _expression.size())
645 0 : break;
646 :
647 : // If we're at the start of a brace delimited section, then we
648 : // parse just that section:
649 0 : if (_expression[nextstart] == '{')
650 : {
651 0 : nextstart++;
652 0 : end = _expression.find('}', nextstart);
653 : }
654 : // otherwise we parse the whole thing
655 : else
656 0 : end = std::string::npos;
657 :
658 : // We either want the whole end of the string (end == npos) or
659 : // a substring in the middle.
660 : _subexpressions.push_back
661 0 : (_expression.substr(nextstart, (end == std::string::npos) ?
662 : std::string::npos : end - nextstart));
663 :
664 : // fparser can crash on empty expressions
665 0 : libmesh_error_msg_if(_subexpressions.back().empty(),
666 : "ERROR: FunctionParser is unable to parse empty expression.\n");
667 :
668 :
669 : #ifdef LIBMESH_HAVE_FPARSER
670 : // Parse (and optimize if possible) the subexpression.
671 : // Add some basic constants, to Real precision.
672 0 : auto fp = std::make_unique<FunctionParserBase<Output>>();
673 0 : fp->AddConstant("NaN", std::numeric_limits<Real>::quiet_NaN());
674 0 : fp->AddConstant("pi", std::acos(Real(-1)));
675 0 : fp->AddConstant("e", std::exp(Real(1)));
676 0 : libmesh_error_msg_if
677 : (fp->Parse(_subexpressions.back(), variables) != -1, // -1 for success
678 : "ERROR: FunctionParser is unable to parse expression: "
679 : << _subexpressions.back() << '\n' << fp->ErrorMsg());
680 0 : fp->Optimize();
681 0 : parsers.push_back(std::move(fp));
682 : #else
683 : libmesh_error_msg("ERROR: This functionality requires fparser!");
684 : #endif
685 :
686 : // If at end, use nextstart=maxSize. Else start at next
687 : // character.
688 0 : nextstart = (end == std::string::npos) ?
689 : std::string::npos : end + 1;
690 : }
691 0 : }
692 :
693 : template <typename Output>
694 : inline
695 : std::size_t
696 0 : ParsedFEMFunction<Output>::find_name (std::string_view varname,
697 : std::string_view expr) const
698 : {
699 0 : const std::size_t namesize = varname.size();
700 0 : std::size_t varname_i = expr.find(varname);
701 :
702 0 : while ((varname_i != std::string::npos) &&
703 0 : (((varname_i > 0) &&
704 0 : (std::isalnum(expr[varname_i-1]) ||
705 0 : (expr[varname_i-1] == '_'))) ||
706 0 : ((varname_i+namesize < expr.size()) &&
707 0 : (std::isalnum(expr[varname_i+namesize]) ||
708 0 : (expr[varname_i+namesize] == '_')))))
709 : {
710 0 : varname_i = expr.find(varname, varname_i+1);
711 : }
712 :
713 0 : return varname_i;
714 : }
715 :
716 :
717 :
718 : // Helper function for evaluating function arguments
719 : template <typename Output>
720 : inline
721 : void
722 0 : ParsedFEMFunction<Output>::eval_args (const FEMContext & c,
723 : const Point & p,
724 : const Real time)
725 : {
726 0 : _spacetime[0] = p(0);
727 : #if LIBMESH_DIM > 1
728 0 : _spacetime[1] = p(1);
729 : #endif
730 : #if LIBMESH_DIM > 2
731 0 : _spacetime[2] = p(2);
732 : #endif
733 0 : _spacetime[LIBMESH_DIM] = time;
734 :
735 0 : unsigned int request_index = 0;
736 0 : for (unsigned int v=0; v != _n_vars; ++v)
737 : {
738 0 : if (!_need_var[v])
739 0 : continue;
740 :
741 0 : c.point_value(v, p, _spacetime[LIBMESH_DIM+1+request_index]);
742 0 : request_index++;
743 : }
744 :
745 0 : if (_n_requested_grad_components)
746 0 : for (unsigned int v=0; v != _n_vars; ++v)
747 : {
748 0 : if (!_need_var_grad[v*LIBMESH_DIM]
749 : #if LIBMESH_DIM > 1
750 0 : && !_need_var_grad[v*LIBMESH_DIM+1]
751 : #if LIBMESH_DIM > 2
752 0 : && !_need_var_grad[v*LIBMESH_DIM+2]
753 : #endif
754 : #endif
755 : )
756 0 : continue;
757 :
758 0 : Gradient g;
759 0 : c.point_gradient(v, p, g);
760 :
761 0 : for (unsigned int d=0; d != LIBMESH_DIM; ++d)
762 : {
763 0 : if (!_need_var_grad[v*LIBMESH_DIM+d])
764 0 : continue;
765 :
766 0 : _spacetime[LIBMESH_DIM+1+request_index] = g(d);
767 0 : request_index++;
768 : }
769 : }
770 :
771 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
772 0 : if (_n_requested_hess_components)
773 0 : for (unsigned int v=0; v != _n_vars; ++v)
774 : {
775 0 : if (!_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM]
776 : #if LIBMESH_DIM > 1
777 0 : && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+1]
778 0 : && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+2]
779 0 : && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+3]
780 : #if LIBMESH_DIM > 2
781 0 : && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+4]
782 0 : && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+5]
783 0 : && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+6]
784 0 : && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+7]
785 0 : && !_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+8]
786 : #endif
787 : #endif
788 : )
789 0 : continue;
790 :
791 0 : Tensor h;
792 0 : c.point_hessian(v, p, h);
793 :
794 0 : for (unsigned int d1=0; d1 != LIBMESH_DIM; ++d1)
795 0 : for (unsigned int d2=0; d2 != LIBMESH_DIM; ++d2)
796 : {
797 0 : if (!_need_var_hess[v*LIBMESH_DIM*LIBMESH_DIM+d1*LIBMESH_DIM+d2])
798 0 : continue;
799 :
800 0 : _spacetime[LIBMESH_DIM+1+request_index] = h(d1,d2);
801 0 : request_index++;
802 : }
803 : }
804 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
805 :
806 0 : if (_requested_normals)
807 : {
808 : FEBase * side_fe;
809 0 : c.get_side_fe(0, side_fe);
810 :
811 0 : const std::vector<Point> & normals = side_fe->get_normals();
812 :
813 0 : const std::vector<Point> & xyz = side_fe->get_xyz();
814 :
815 0 : libmesh_assert_equal_to(normals.size(), xyz.size());
816 :
817 : // We currently only support normals at quadrature points!
818 : #ifndef NDEBUG
819 0 : bool at_quadrature_point = false;
820 : #endif
821 0 : for (auto qp : index_range(normals))
822 : {
823 0 : if (p == xyz[qp])
824 : {
825 0 : const Point & n = normals[qp];
826 0 : for (unsigned int d=0; d != LIBMESH_DIM; ++d)
827 : {
828 0 : _spacetime[LIBMESH_DIM+1+request_index] = n(d);
829 0 : request_index++;
830 : }
831 : #ifndef NDEBUG
832 0 : at_quadrature_point = true;
833 : #endif
834 0 : break;
835 : }
836 : }
837 :
838 0 : libmesh_assert(at_quadrature_point);
839 : }
840 :
841 : // The remaining locations in _spacetime are currently set at construction
842 : // but could potentially be made dynamic
843 0 : }
844 :
845 :
846 : // Evaluate the ith FunctionParser and check the result
847 : #ifdef LIBMESH_HAVE_FPARSER
848 : template <typename Output>
849 : inline
850 : Output
851 0 : ParsedFEMFunction<Output>::eval (FunctionParserBase<Output> & parser,
852 : std::string_view libmesh_dbg_var(function_name),
853 : unsigned int libmesh_dbg_var(component_idx)) const
854 : {
855 : #ifndef NDEBUG
856 0 : Output result = parser.Eval(_spacetime.data());
857 0 : int error_code = parser.EvalError();
858 0 : if (error_code)
859 : {
860 0 : libMesh::err << "ERROR: FunctionParser is unable to evaluate component "
861 0 : << component_idx
862 0 : << " for '"
863 0 : << function_name;
864 :
865 0 : for (auto i : index_range(parsers))
866 0 : if (parsers[i].get() == &parser)
867 0 : libMesh::err << "' of expression '"
868 0 : << _subexpressions[i];
869 :
870 0 : libMesh::err << "' with arguments:\n";
871 0 : for (const auto & st : _spacetime)
872 0 : libMesh::err << '\t' << st << '\n';
873 0 : libMesh::err << '\n';
874 :
875 : // Currently no API to report error messages, we'll do it manually
876 0 : std::string error_message = "Reason: ";
877 :
878 0 : switch (error_code)
879 : {
880 0 : case 1:
881 0 : error_message += "Division by zero";
882 0 : break;
883 0 : case 2:
884 0 : error_message += "Square Root error (negative value)";
885 0 : break;
886 0 : case 3:
887 0 : error_message += "Log error (negative value)";
888 0 : break;
889 0 : case 4:
890 0 : error_message += "Trigonometric error (asin or acos of illegal value)";
891 0 : break;
892 0 : case 5:
893 0 : error_message += "Maximum recursion level reached";
894 0 : break;
895 0 : default:
896 0 : error_message += "Unknown";
897 0 : break;
898 : }
899 0 : libmesh_error_msg(error_message);
900 : }
901 :
902 0 : return result;
903 : #else
904 0 : return parser.Eval(_spacetime.data());
905 : #endif
906 : }
907 : #else // LIBMESH_HAVE_FPARSER
908 : template <typename Output>
909 : inline
910 : Output
911 : ParsedFEMFunction<Output>::eval (char & /*parser*/,
912 : std::string_view /*function_name*/,
913 : unsigned int /*component_idx*/) const
914 : {
915 : libmesh_error_msg("ERROR: This functionality requires fparser!");
916 : return Output(0);
917 : }
918 : #endif
919 :
920 :
921 : } // namespace libMesh
922 :
923 : #endif // LIBMESH_PARSED_FEM_FUNCTION_H
|
