Reference documentation for deal.II version Git f40be01994 2020-04-09 07:13:12 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
symengine_number_types.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 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 
20 # include <deal.II/base/signaling_nan.h>
21 
22 # include <deal.II/differentiation/sd/symengine_number_types.h>
23 # include <deal.II/differentiation/sd/symengine_types.h>
24 # include <deal.II/differentiation/sd/symengine_utilities.h>
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 
36 DEAL_II_NAMESPACE_OPEN
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,
149  Utilities::convert_expression_vector_to_basic_vector(arguments)))
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 types::substitution_map &substitution_values) const
288  {
289  return Expression(get_expression().subs(
290  Utilities::convert_expression_map_to_basic_map(substitution_values)));
291  }
292 
293 
294  Expression
296  const Expression &value) const
297  {
298  Assert(SE::is_a<SE::Symbol>(symbol.get_value()),
299  ExcMessage(
300  "Substitution with a number that does not represent a symbol."));
301 
302  types::substitution_map sub_vals;
303  sub_vals[symbol] = value;
304  return substitute(sub_vals);
305  }
306 
307 
308  /* -------------------- Math and relational operators ------------------ */
309 
310 
311  Expression &
313  {
314  if (this != &rhs)
315  this->expression = rhs.get_expression();
316 
317  return *this;
318  }
319 
320 
321  Expression &
323  {
324  if (this != &rhs)
325  this->expression = std::move(rhs.expression);
326 
327  return *this;
328  }
329 
330 
331  Expression
333  {
334  return Expression(-get_expression());
335  }
336 
337 
338  Expression &
340  {
341  this->expression += rhs.get_expression();
342  return *this;
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  std::ostream &
371  operator<<(std::ostream &stream, const Expression &expr)
372  {
373  stream << expr.get_expression();
374  return stream;
375  }
376 
377 
378  std::istream &
379  operator>>(std::istream &stream, Expression &expr)
380  {
381  std::string str;
382  stream >> str;
383  expr.parse(str);
384  return stream;
385  }
386 
387 
388  Expression
389  operator==(const Expression &lhs, const Expression &rhs)
390  {
391  return Expression(SE::Eq(lhs.get_RCP(), rhs.get_RCP()));
392  }
393 
394 
395  Expression
396  operator!=(const Expression &lhs, const Expression &rhs)
397  {
398  return Expression(SE::Ne(lhs.get_RCP(), rhs.get_RCP()));
399  }
400 
401 
402  Expression
403  operator<(const Expression &lhs, const Expression &rhs)
404  {
405  return Expression(SE::Lt(lhs.get_RCP(), rhs.get_RCP()));
406  }
407 
408 
409  Expression
410  operator>(const Expression &lhs, const Expression &rhs)
411  {
412  return Expression(SE::Gt(lhs.get_RCP(), rhs.get_RCP()));
413  }
414 
415 
416  Expression
417  operator<=(const Expression &lhs, const Expression &rhs)
418  {
419  return Expression(SE::Le(lhs.get_RCP(), rhs.get_RCP()));
420  }
421 
422 
423  Expression
424  operator>=(const Expression &lhs, const Expression &rhs)
425  {
426  return Expression(SE::Ge(lhs.get_RCP(), rhs.get_RCP()));
427  }
428 
429 
431  {
432  Assert(SE::is_a_Boolean(expression.get_value()),
433  ExcMessage("The expression must return a boolean type."));
434 
435  const SE::RCP<const SE::Boolean> expression_rcp =
436  SE::rcp_static_cast<const SE::Boolean>(expression.get_RCP());
437 
438  return Expression(SE::logical_not(expression_rcp));
439  }
440 
441 
442  Expression operator&(const Expression &lhs, const Expression &rhs)
443  {
444  Assert(SE::is_a_Boolean(lhs.get_value()),
445  ExcMessage("The lhs expression must return a boolean type."));
446  Assert(SE::is_a_Boolean(rhs.get_value()),
447  ExcMessage("The rhs expression must return a boolean type."));
448 
449  const SE::RCP<const SE::Boolean> lhs_rcp =
450  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
451  const SE::RCP<const SE::Boolean> rhs_rcp =
452  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
453 
454  return Expression(SE::logical_and({lhs_rcp, rhs_rcp}));
455  }
456 
457 
458  Expression
459  operator|(const Expression &lhs, const Expression &rhs)
460  {
461  Assert(SE::is_a_Boolean(lhs.get_value()),
462  ExcMessage("The lhs expression must return a boolean type."));
463  Assert(SE::is_a_Boolean(rhs.get_value()),
464  ExcMessage("The rhs expression must return a boolean type."));
465 
466  const SE::RCP<const SE::Boolean> lhs_rcp =
467  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
468  const SE::RCP<const SE::Boolean> rhs_rcp =
469  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
470 
471  return Expression(SE::logical_or({lhs_rcp, rhs_rcp}));
472  }
473 
474 
475  Expression
476  operator^(const Expression &lhs, const Expression &rhs)
477  {
478  Assert(SE::is_a_Boolean(lhs.get_value()),
479  ExcMessage("The lhs expression must return a boolean type."));
480  Assert(SE::is_a_Boolean(rhs.get_value()),
481  ExcMessage("The rhs expression must return a boolean type."));
482 
483  const SE::RCP<const SE::Boolean> lhs_rcp =
484  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
485  const SE::RCP<const SE::Boolean> rhs_rcp =
486  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
487 
488  return Expression(SE::logical_xor({lhs_rcp, rhs_rcp}));
489  }
490 
491 
492  Expression
493  operator&&(const Expression &lhs, const Expression &rhs)
494  {
495  return lhs & rhs;
496  }
497 
498 
499  Expression
500  operator||(const Expression &lhs, const Expression &rhs)
501  {
502  return lhs | rhs;
503  }
504 
505 
506  Expression
507  operator+(Expression lhs, const Expression &rhs)
508  {
509  lhs += rhs;
510  return lhs;
511  }
512 
513 
514  Expression
515  operator-(Expression lhs, const Expression &rhs)
516  {
517  lhs -= rhs;
518  return lhs;
519  }
520 
521 
523  {
524  lhs *= rhs;
525  return lhs;
526  }
527 
528 
529  Expression
530  operator/(Expression lhs, const Expression &rhs)
531  {
532  lhs /= rhs;
533  return lhs;
534  }
535 
536 
537  } // namespace SD
538 } // namespace Differentiation
539 
540 
541 DEAL_II_NAMESPACE_CLOSE
542 
543 #endif // DEAL_II_WITH_SYMENGINE
std::ostream & operator<<(std::ostream &stream, const Expression &expression)
Expression operator+(Expression lhs, const Expression &rhs)
Expression operator<(const Expression &lhs, const Expression &rhs)
Expression differentiate(const Expression &symbol) const
std::istream & operator>>(std::istream &stream, Expression &expression)
Expression & operator/=(const Expression &rhs)
Expression operator &&(const Expression &lhs, const Expression &rhs)
Expression operator>=(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcIO()
Expression & parse(const std::string &expression)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
std::ostream & print(std::ostream &stream) const
Expression substitute(const types::substitution_map &substitution_values) const
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
Expression operator/(Expression lhs, const Expression &rhs)
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression operator<=(const Expression &lhs, const Expression &rhs)
#define Assert(cond, exc)
Definition: exceptions.h:1419
Expression & operator*=(const Expression &rhs)
Expression operator|(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcSymEngineParserError(std::string arg1)
const SymEngine::Basic & get_value() const
const SymEngine::Expression & get_expression() const
Expression & operator=(const Expression &rhs)
Expression operator>(const Expression &lhs, const Expression &rhs)
void save(std::ostream &stream) const
Expression operator*(Expression lhs, const Expression &rhs)
Expression operator!(const Expression &expression)
Definition: cuda.h:31
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 &rhs)
Expression & operator+=(const Expression &rhs)
Expression operator!=(const Expression &lhs, const Expression &rhs)
Expression operator &(const Expression &lhs, const Expression &rhs)