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