Reference documentation for deal.II version GIT 6a72d26406 2023-06-07 13:05:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
symengine_number_types.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/config.h>
17 
18 #ifdef DEAL_II_WITH_SYMENGINE
19 
21 
25 
26 # include <symengine/complex_double.h>
27 # include <symengine/logic.h>
28 # include <symengine/number.h>
29 # include <symengine/parser.h>
30 # include <symengine/real_double.h>
31 # include <symengine/symbol.h>
32 # include <symengine/symengine_exception.h>
33 
34 # include <string>
35 
37 
38 namespace Differentiation
39 {
40  namespace SD
41  {
42  namespace SE = ::SymEngine;
43 
44 
45  /* ---------------------------- Constructors -------------------------- */
46 
47 
49  : expression()
50  {}
51 
52 
53  Expression::Expression(const bool value)
54  : expression(SE::boolean(value))
55  {}
56 
57 
58  Expression::Expression(const SymEngine::integer_class &value)
59  : expression(value)
60  {}
61 
62 
63  Expression::Expression(const SymEngine::rational_class &value)
64  : expression(value)
65  {}
66 
67 
69  const Expression &expression_if_true,
70  const Expression &expression_if_false)
71  {
72  Assert(SE::is_a_Boolean(condition.get_value()),
73  ExcMessage(
74  "The conditional expression must return a boolean type."));
75 
76  const SE::RCP<const SE::Boolean> condition_rcp =
77  SE::rcp_static_cast<const SE::Boolean>(condition.get_RCP());
78  expression =
79  SE::piecewise({{expression_if_true.get_RCP(), condition_rcp},
80  {expression_if_false.get_RCP(), SE::boolTrue}});
81  }
82 
83 
84  Expression::Expression(const std::vector<std::pair<Expression, Expression>>
85  & condition_expression,
86  const Expression &expression_otherwise)
87  {
88  SE::PiecewiseVec piecewise_function;
89  piecewise_function.reserve(condition_expression.size() + 1);
90 
91  // Add tested conditional entries
92  for (const auto &entry : condition_expression)
93  {
94  Assert(SE::is_a_Boolean(entry.first.get_value()),
95  ExcMessage(
96  "The conditional expression must return a boolean type."));
97  piecewise_function.push_back(
98  {entry.second.get_RCP(),
99  SE::rcp_static_cast<const SE::Boolean>(entry.first.get_RCP())});
100  }
101 
102  // Add default value
103  piecewise_function.push_back(
104  {expression_otherwise.get_RCP(), SE::boolTrue});
105 
106  // Initialize
107  expression = SE::piecewise(std::move(piecewise_function));
108  }
109 
110 
111  Expression::Expression(const std::vector<std::pair<Expression, Expression>>
112  &condition_expression)
113  {
114  // Use the other constructor with a fatal termination point
115  // ensuring that an error is thrown if none of the conditions
116  // are met.
117  *this = Expression(condition_expression,
118  Expression(numbers::signaling_nan<double>()));
119  }
120 
121 
122  Expression::Expression(const char *symbol)
123  : expression(SE::symbol(symbol))
124  {}
125 
126 
127  Expression::Expression(const std::string &str,
128  const bool parse_as_expression)
129  {
130  try
131  {
132  expression = (parse_as_expression ?
133  SE::parse(str) // The string is a symbolic "name"
134  :
135  SE::rcp_static_cast<const SE::Basic>(SE::symbol(
136  str))); // The string is a symbolic "expression"
137  }
138  catch (...)
139  {
141  }
142  }
143 
144 
145  Expression::Expression(const std::string & symbol_func,
146  const types::symbol_vector &arguments)
147  : expression(SE::function_symbol(
148  symbol_func,
150  {}
151 
152 
153  Expression::Expression(const SymEngine::Expression &rhs)
154  : expression(rhs)
155  {}
156 
157 
158  Expression::Expression(const SymEngine::RCP<const SymEngine::Basic> &rhs)
159  : expression(rhs)
160  {}
161 
162 
163  Expression::Expression(SymEngine::RCP<const SymEngine::Basic> &&rhs)
164  : expression(rhs)
165  {}
166 
167 
168  /* ------------------------------ Utilities ---------------------------- */
169 
170 
171  Expression &
172  Expression::parse(const std::string &expression)
173  {
174  *this = Expression(expression, true /*parse_as_expression*/);
175  return *this;
176  }
177 
178 
179  std::ostream &
180  Expression::print(std::ostream &os) const
181  {
182  os << *this;
183  return os;
184  }
185 
186 
187  void
188  Expression::save(std::ostream &os) const
189  {
190  // We write each expression on a new line.
191  // Note: SymEngine outputs a non-terminating string
192  os << *this;
193  os << std::endl;
194  }
195 
196 
197  void
198  Expression::load(std::istream &is)
199  {
200  // Need to make sure that we read the entire line in,
201  // and then subsequently parse it.
202  std::string expr;
203  std::getline(is, expr);
204  Assert(!is.bad(), ExcIO());
205  parse(expr);
206  }
207 
208 
209  /* ------------------------------- Values ----------------------------- */
210 
211 
212  const SE::Expression &
214  {
215  return expression;
216  }
217 
218 
219  SE::Expression &
221  {
222  return expression;
223  }
224 
225 
226  const SE::Basic &
228  {
229  return *get_RCP();
230  }
231 
232 
233  const SE::RCP<const SE::Basic> &
235  {
236  return get_expression().get_basic();
237  }
238 
239 
240  /* --------------------------- Differentiation ------------------------- */
241 
242 
243  Expression
245  const SymEngine::RCP<const SymEngine::Symbol> &symbol) const
246  {
247  return Expression(SE::diff(get_RCP(), symbol));
248  }
249 
250 
251  Expression
253  const SymEngine::RCP<const SymEngine::Basic> &symbol) const
254  {
255  // Potential symbol
256  return Expression(SE::sdiff(get_RCP(), symbol));
257  }
258 
259 
260  Expression
262  {
263  return differentiate(symbol.get_RCP());
264  }
265 
266 
267  /* ------------- Conversion operators ------------------------- */
268 
269 
270  Expression::operator const SymEngine::Expression &() const
271  {
272  return get_expression();
273  }
274 
275 
276  Expression::operator const SymEngine::RCP<const SymEngine::Basic> &() const
277  {
278  return get_expression().get_basic();
279  }
280 
281 
282  /* ------------- Dictionary-based substitution ------------------------- */
283 
284 
285  Expression
287  const SymEngine::map_basic_basic &substitution_values) const
288  {
289  return Expression(get_expression().subs(substitution_values));
290  }
291 
292 
293  Expression
295  const types::substitution_map &substitution_values) const
296  {
297  return substitute(
299  }
300 
301 
302  Expression
304  const Expression &value) const
305  {
306  Assert(SE::is_a<SE::Symbol>(symbol.get_value()),
307  ExcMessage(
308  "Substitution with a number that does not represent a symbol."));
309 
310  types::substitution_map sub_vals;
311  sub_vals[symbol] = value;
312  return substitute(sub_vals);
313  }
314 
315 
316  /* -------------------- Math and relational operators ------------------ */
317 
318 
319  Expression &
321  {
322  if (this != &rhs)
323  this->expression = rhs.get_expression();
324 
325  return *this;
326  }
327 
328 
329  Expression &
331  {
332  if (this != &rhs)
333  this->expression = std::move(rhs.expression);
334 
335  return *this;
336  }
337 
338 
339  Expression
341  {
342  return Expression(-get_expression());
343  }
344 
345 
346  Expression &
348  {
349  this->expression += rhs.get_expression();
350  return *this;
351  }
352 
353 
354  Expression &
356  {
357  this->expression -= rhs.get_expression();
358  return *this;
359  }
360 
361 
362  Expression &
364  {
365  this->expression *= rhs.get_expression();
366  return *this;
367  }
368 
369 
370  Expression &
372  {
373  this->expression /= rhs.get_expression();
374  return *this;
375  }
376 
377 
378  std::ostream &
379  operator<<(std::ostream &stream, const Expression &expr)
380  {
381  stream << expr.get_expression();
382  return stream;
383  }
384 
385 
386  std::istream &
387  operator>>(std::istream &stream, Expression &expr)
388  {
389  std::string str;
390  stream >> str;
391  expr.parse(str);
392  return stream;
393  }
394 
395 
396  Expression
397  operator==(const Expression &lhs, const Expression &rhs)
398  {
399  return Expression(SE::Eq(lhs.get_RCP(), rhs.get_RCP()));
400  }
401 
402 
403  Expression
404  operator!=(const Expression &lhs, const Expression &rhs)
405  {
406  return Expression(SE::Ne(lhs.get_RCP(), rhs.get_RCP()));
407  }
408 
409 
411  operator<(const Expression &lhs, const Expression &rhs)
412  {
413  return Expression(SE::Lt(lhs.get_RCP(), rhs.get_RCP()));
414  }
415 
416 
417  Expression
418  operator>(const Expression &lhs, const Expression &rhs)
419  {
420  return Expression(SE::Gt(lhs.get_RCP(), rhs.get_RCP()));
421  }
422 
423 
425  operator<=(const Expression &lhs, const Expression &rhs)
426  {
427  return Expression(SE::Le(lhs.get_RCP(), rhs.get_RCP()));
428  }
429 
430 
431  Expression
432  operator>=(const Expression &lhs, const Expression &rhs)
433  {
434  return Expression(SE::Ge(lhs.get_RCP(), rhs.get_RCP()));
435  }
436 
437 
438  Expression
439  operator!(const Expression &expression)
440  {
441  Assert(SE::is_a_Boolean(expression.get_value()),
442  ExcMessage("The expression must return a boolean type."));
443 
444  const SE::RCP<const SE::Boolean> expression_rcp =
445  SE::rcp_static_cast<const SE::Boolean>(expression.get_RCP());
446 
447  return Expression(SE::logical_not(expression_rcp));
448  }
449 
450 
451  Expression
452  operator&(const Expression &lhs, const Expression &rhs)
453  {
454  Assert(SE::is_a_Boolean(lhs.get_value()),
455  ExcMessage("The lhs expression must return a boolean type."));
456  Assert(SE::is_a_Boolean(rhs.get_value()),
457  ExcMessage("The rhs expression must return a boolean type."));
458 
459  const SE::RCP<const SE::Boolean> lhs_rcp =
460  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
461  const SE::RCP<const SE::Boolean> rhs_rcp =
462  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
463 
464  return Expression(SE::logical_and({lhs_rcp, rhs_rcp}));
465  }
466 
467 
468  Expression
469  operator|(const Expression &lhs, const Expression &rhs)
470  {
471  Assert(SE::is_a_Boolean(lhs.get_value()),
472  ExcMessage("The lhs expression must return a boolean type."));
473  Assert(SE::is_a_Boolean(rhs.get_value()),
474  ExcMessage("The rhs expression must return a boolean type."));
475 
476  const SE::RCP<const SE::Boolean> lhs_rcp =
477  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
478  const SE::RCP<const SE::Boolean> rhs_rcp =
479  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
480 
481  return Expression(SE::logical_or({lhs_rcp, rhs_rcp}));
482  }
483 
484 
485  Expression
486  operator^(const Expression &lhs, const Expression &rhs)
487  {
488  Assert(SE::is_a_Boolean(lhs.get_value()),
489  ExcMessage("The lhs expression must return a boolean type."));
490  Assert(SE::is_a_Boolean(rhs.get_value()),
491  ExcMessage("The rhs expression must return a boolean type."));
492 
493  const SE::RCP<const SE::Boolean> lhs_rcp =
494  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
495  const SE::RCP<const SE::Boolean> rhs_rcp =
496  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
497 
498  return Expression(SE::logical_xor({lhs_rcp, rhs_rcp}));
499  }
500 
501 
502  Expression
503  operator&&(const Expression &lhs, const Expression &rhs)
504  {
505  return lhs & rhs;
506  }
507 
508 
509  Expression
510  operator||(const Expression &lhs, const Expression &rhs)
511  {
512  return lhs | rhs;
513  }
514 
515 
516  Expression
517  operator+(Expression lhs, const Expression &rhs)
518  {
519  lhs += rhs;
520  return lhs;
521  }
522 
523 
524  Expression
525  operator-(Expression lhs, const Expression &rhs)
526  {
527  lhs -= rhs;
528  return lhs;
529  }
530 
531 
532  Expression
533  operator*(Expression lhs, const Expression &rhs)
534  {
535  lhs *= rhs;
536  return lhs;
537  }
538 
539 
540  Expression
541  operator/(Expression lhs, const Expression &rhs)
542  {
543  lhs /= rhs;
544  return lhs;
545  }
546 
547 
548  } // namespace SD
549 } // namespace Differentiation
550 
551 
553 
554 #endif // DEAL_II_WITH_SYMENGINE
Expression & operator/=(const Expression &rhs)
Expression & parse(const std::string &expression)
Expression & operator=(const Expression &rhs)
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
Expression & operator*=(const Expression &rhs)
std::ostream & print(std::ostream &stream) const
Expression substitute(const types::substitution_map &substitution_values) const
void save(std::ostream &stream) const
Expression & operator-=(const Expression &rhs)
const SymEngine::Basic & get_value() const
const SymEngine::Expression & get_expression() const
Expression & operator+=(const Expression &rhs)
Expression differentiate(const Expression &symbol) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcSymEngineParserError(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
SymEngine::map_basic_basic convert_expression_map_to_basic_map(const SD::types::substitution_map &substitution_map)
SymEngine::vec_basic convert_expression_vector_to_basic_vector(const SD::types::symbol_vector &symbol_vector)
std::vector< SD::Expression > symbol_vector
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
Expression operator*(Expression lhs, const Expression &rhs)
Expression operator<(const Expression &lhs, const Expression &rhs)
Expression operator-(Expression lhs, const Expression &rhs)
Expression operator+(Expression lhs, const Expression &rhs)
Expression operator!(const Expression &expression)
Expression operator||(const Expression &lhs, const Expression &rhs)
Expression operator^(const Expression &lhs, const Expression &rhs)
Expression operator>=(const Expression &lhs, const Expression &rhs)
Expression operator!=(const Expression &lhs, const Expression &rhs)
Expression operator|(const Expression &lhs, const Expression &rhs)
Expression operator>(const Expression &lhs, const Expression &rhs)
Expression operator&(const Expression &lhs, const Expression &rhs)
Expression operator/(Expression lhs, const Expression &rhs)
Expression operator<=(const Expression &lhs, const Expression &rhs)
Expression operator&&(const Expression &lhs, const Expression &rhs)
std::ostream & operator<<(std::ostream &stream, const Expression &expression)
Expression operator==(const Expression &lhs, const Expression &rhs)
std::istream & operator>>(std::istream &stream, Expression &expression)
T logical_or(const T &t, const MPI_Comm mpi_communicator)