Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
38namespace 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()),
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());
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()),
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
245 const SymEngine::RCP<const SymEngine::Symbol> &symbol) const
246 {
247 return Expression(SE::diff(get_RCP(), symbol));
248 }
249
250
253 const SymEngine::RCP<const SymEngine::Basic> &symbol) const
254 {
255 // Potential symbol
256 return Expression(SE::sdiff(get_RCP(), symbol));
257 }
258
259
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
295 const types::substitution_map &substitution_values) const
296 {
297 return substitute(
299 }
300
301
304 const Expression &value) const
305 {
306 Assert(SE::is_a<SE::Symbol>(symbol.get_value()),
308 "Substitution with a number that does not represent a symbol."));
309
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
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
518 {
519 lhs += rhs;
520 return lhs;
521 }
522
523
524 Expression
526 {
527 lhs -= rhs;
528 return lhs;
529 }
530
531
532 Expression
534 {
535 lhs *= rhs;
536 return lhs;
537 }
538
539
540 Expression
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcSymEngineParserError(std::string arg1)
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
SymEngine::map_basic_basic convert_expression_map_to_basic_map(const SD::types::substitution_map &substitution_map)
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)