Reference documentation for deal.II version 9.4.1
\(\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.h
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#ifndef dealii_differentiation_sd_symengine_number_types_h
17#define dealii_differentiation_sd_symengine_number_types_h
18
19#include <deal.II/base/config.h>
20
21#ifdef DEAL_II_WITH_SYMENGINE
22
23// Low level
24# include <symengine/basic.h>
25# include <symengine/dict.h>
26# include <symengine/symengine_exception.h>
27# include <symengine/symengine_rcp.h>
28
29// Number types
30# include <symengine/expression.h>
31# include <symengine/integer.h>
32# include <symengine/logic.h>
33# include <symengine/number.h>
34# include <symengine/rational.h>
35
36// Number operations
37# include <symengine/add.h>
38# include <symengine/functions.h>
39# include <symengine/mul.h>
40# include <symengine/pow.h>
41
42// Evaluation
43# include <symengine/eval.h>
44# include <symengine/eval_arb.h>
45# include <symengine/eval_double.h>
46
47// Differentiation
49# include <deal.II/base/numbers.h>
50
53
54# include <boost/serialization/split_member.hpp>
55
56# include <symengine/derivative.h>
57
58# include <algorithm>
59# include <memory>
60# include <sstream>
61# include <type_traits>
62# include <utility>
63# include <vector>
64
66
67namespace Differentiation
68{
69 namespace SD
70 {
82 std::string,
83 << "The string '" << arg1
84 << "' could not be parsed successfully. Are you sure that (1) it "
85 << "consists of legitimate operations and syntax, and (2) you've "
86 << "previously declared all symbolic variables that are present "
87 << "in the expression?");
88
90
178 {
179 public:
184
188 Expression();
189
197 explicit Expression(const bool value);
198
206 template <typename NumberType,
207 typename = typename std::enable_if<
208 std::is_arithmetic<NumberType>::value>::type>
209 explicit Expression(const NumberType &value);
210
218 template <typename NumberType,
219 typename = typename std::enable_if<
220 std::is_arithmetic<NumberType>::value>::type>
221 explicit Expression(const std::complex<NumberType> &value);
222
226 Expression(const SymEngine::integer_class &value);
227
234 template <typename NumberType,
235 typename = typename std::enable_if<
236 std::is_integral<NumberType>::value>::type>
237 Expression(const NumberType &numerator, const NumberType &denominator);
238
242 Expression(const SymEngine::rational_class &value);
243
280 Expression(const Expression &condition,
281 const Expression &expression_if_true,
282 const Expression &expression_if_false);
283
305 Expression(const std::vector<std::pair<Expression, Expression>>
306 & condition_expression,
307 const Expression &expression_otherwise);
308
328 Expression(const std::vector<std::pair<Expression, Expression>>
329 &condition_expression);
330
331
338 Expression(const char *symbol);
339
352 Expression(const std::string &symb_expr,
353 const bool parse_as_expression = false);
354
361 Expression(const std::string & symbol_func,
362 const types::symbol_vector &arguments);
363
367 Expression(const Expression &rhs) = default;
368
376 explicit Expression(const SymEngine::Expression &rhs);
377
386 Expression(const SymEngine::RCP<const SymEngine::Basic> &rhs);
387
391 Expression(Expression &&rhs) = default;
392
401 Expression(SymEngine::RCP<const SymEngine::Basic> &&rhs);
402
406 virtual ~Expression() = default;
407
409
414
420 Expression &
421 parse(const std::string &expression);
422
429 std::ostream &
430 print(std::ostream &stream) const;
431
437 void
438 save(std::ostream &stream) const;
439
451 void
452 load(std::istream &stream);
453
462 template <class Archive>
463 void
464 save(Archive &archive, const unsigned int version) const;
465
478 template <class Archive>
479 void
480 load(Archive &archive, const unsigned int version);
481
482# ifdef DOXYGEN
497 template <class Archive>
498 void
499 serialize(Archive &archive, const unsigned int version);
500# else
501 // This macro defines the serialize() method that is compatible with
502 // the templated save() and load() method that have been implemented.
503 BOOST_SERIALIZATION_SPLIT_MEMBER()
504# endif
505
507
512
516 const SymEngine::Expression &
517 get_expression() const;
518
523 const SymEngine::Basic &
524 get_value() const;
525
530 const SymEngine::RCP<const SymEngine::Basic> &
531 get_RCP() const;
532
534
539
546 Expression &
547 operator=(const Expression &rhs);
548
555 Expression &
556 operator=(Expression &&rhs) noexcept;
557
564 Expression &
565 operator+=(const Expression &rhs);
566
573 Expression &
574 operator-=(const Expression &rhs);
575
582 Expression &
583 operator*=(const Expression &rhs);
584
591 Expression &
592 operator/=(const Expression &rhs);
593
595
600
607 template <typename NumberType>
608 Expression &
609 operator=(const NumberType &rhs);
610
620 operator-() const;
621
628 template <typename NumberType>
629 Expression &
630 operator+=(const NumberType &rhs);
631
638 template <typename NumberType>
639 Expression &
640 operator-=(const NumberType &rhs);
641
648 template <typename NumberType>
649 Expression &
650 operator*=(const NumberType &rhs);
651
658 template <typename NumberType>
659 Expression &
660 operator/=(const NumberType &rhs);
661
663
668
674 differentiate(const Expression &symbol) const;
675
682 const SymEngine::RCP<const SymEngine::Symbol> &symbol) const;
683
689 differentiate(const SymEngine::RCP<const SymEngine::Basic> &symbol) const;
690
692
697
711 substitute(const types::substitution_map &substitution_values) const;
712
730 substitute(const SymEngine::map_basic_basic &substitution_values) const;
731
742 substitute(const Expression &symbol, const Expression &value) const;
743
753 template <typename NumberType>
755 substitute(const Expression &symbol, const NumberType &value) const;
756
764 template <typename ReturnType>
765 ReturnType
767 const types::substitution_map &substitution_values) const;
768
780 template <typename ReturnType>
781 ReturnType
783 const SymEngine::map_basic_basic &substitution_values) const;
784
786
791
842 template <typename ResultType>
843 explicit operator ResultType() const;
844
849 explicit operator const SymEngine::Expression &() const;
850
855 operator const SymEngine::RCP<const SymEngine::Basic> &() const;
856
858
859 protected:
863 SymEngine::Expression &
865
866 private:
871 SymEngine::Expression expression;
872 };
873
878
884 template <>
885 struct is_symengine_number<Expression> : std::true_type
886 {};
887
888
894 template <>
895 struct is_sd_number<Expression> : std::true_type
896 {};
897
899
904
910 std::ostream &
911 operator<<(std::ostream &stream, const Expression &expression);
912
918 std::istream &
919 operator>>(std::istream &stream, Expression &expression);
920
922
927
934 operator==(const Expression &lhs, const Expression &rhs);
935
942 operator!=(const Expression &lhs, const Expression &rhs);
943
950 operator<(const Expression &lhs, const Expression &rhs);
951
958 operator>(const Expression &lhs, const Expression &rhs);
959
966 operator<=(const Expression &lhs, const Expression &rhs);
967
974 operator>=(const Expression &lhs, const Expression &rhs);
975
977
982
990 operator!(const Expression &expression);
991
999 operator&(const Expression &lhs, const Expression &rhs);
1000
1008 operator|(const Expression &lhs, const Expression &rhs);
1009
1017 operator^(const Expression &lhs, const Expression &rhs);
1018
1027 operator&&(const Expression &lhs, const Expression &rhs);
1028
1037 operator||(const Expression &lhs, const Expression &rhs);
1038
1040
1045
1052 operator+(Expression lhs, const Expression &rhs);
1053
1060 operator-(Expression lhs, const Expression &rhs);
1061
1068 operator*(Expression lhs, const Expression &rhs);
1069
1076 operator/(Expression lhs, const Expression &rhs);
1077
1086 template <typename NumberType,
1087 typename = typename std::enable_if<
1088 std::is_constructible<Expression, NumberType>::value>::type>
1089 inline Expression
1090 operator+(const NumberType &lhs, const Expression &rhs)
1091 {
1092 return Expression(lhs) + rhs;
1093 }
1094
1103 template <typename NumberType,
1104 typename = typename std::enable_if<
1105 std::is_constructible<Expression, NumberType>::value>::type>
1106 inline Expression
1107 operator+(const Expression &lhs, const NumberType &rhs)
1108 {
1109 return lhs + Expression(rhs);
1110 }
1111
1120 template <typename NumberType,
1121 typename = typename std::enable_if<
1122 std::is_constructible<Expression, NumberType>::value>::type>
1123 inline Expression
1124 operator-(const NumberType &lhs, const Expression &rhs)
1125 {
1126 return Expression(lhs) - rhs;
1127 }
1128
1137 template <typename NumberType,
1138 typename = typename std::enable_if<
1139 std::is_constructible<Expression, NumberType>::value>::type>
1140 inline Expression
1141 operator-(const Expression &lhs, const NumberType &rhs)
1142 {
1143 return lhs - Expression(rhs);
1144 }
1145
1154 template <typename NumberType,
1155 typename = typename std::enable_if<
1156 std::is_constructible<Expression, NumberType>::value>::type>
1157 inline Expression
1158 operator*(const NumberType &lhs, const Expression &rhs)
1159 {
1160 return Expression(lhs) * rhs;
1161 }
1162
1171 template <typename NumberType,
1172 typename = typename std::enable_if<
1173 std::is_constructible<Expression, NumberType>::value>::type>
1174 inline Expression
1175 operator*(const Expression &lhs, const NumberType &rhs)
1176 {
1177 return lhs * Expression(rhs);
1178 }
1179
1188 template <typename NumberType,
1189 typename = typename std::enable_if<
1190 std::is_constructible<Expression, NumberType>::value>::type>
1191 inline Expression
1192 operator/(const NumberType &lhs, const Expression &rhs)
1193 {
1194 return Expression(lhs) / rhs;
1195 }
1196
1205 template <typename NumberType,
1206 typename = typename std::enable_if<
1207 std::is_constructible<Expression, NumberType>::value>::type>
1208 inline Expression
1209 operator/(const Expression &lhs, const NumberType &rhs)
1210 {
1211 return lhs / Expression(rhs);
1212 }
1213
1215
1216 } // namespace SD
1217} // namespace Differentiation
1218
1219
1220/* ---------------------- inline and template functions -------------------- */
1221
1222
1223# ifndef DOXYGEN
1224
1225
1226namespace Differentiation
1227{
1228 namespace SD
1229 {
1230 template <typename NumberType, typename>
1231 Expression::Expression(const NumberType &value)
1232 : expression(value)
1233 {}
1234
1235
1236 template <typename NumberType, typename>
1237 Expression::Expression(const std::complex<NumberType> &value)
1238 : expression(value)
1239 {}
1240
1241
1242 template <typename NumberType, typename>
1243 Expression::Expression(const NumberType &numerator,
1244 const NumberType &denominator)
1245 : expression(
1246 SymEngine::Rational::from_two_ints(*SymEngine::integer(numerator),
1247 *SymEngine::integer(denominator)))
1248 {}
1249
1250
1251 template <class Archive>
1252 void
1253 Expression::save(Archive &ar, const unsigned int /*version*/) const
1254 {
1255 std::stringstream sstream;
1256 sstream << *this;
1257 const std::string expr = sstream.str();
1258 ar & expr;
1259 }
1260
1261
1262 template <class Archive>
1263 void
1264 Expression::load(Archive &ar, const unsigned int /*version*/)
1265 {
1266 std::string expr;
1267 ar & expr;
1268 parse(expr);
1269 }
1270
1271
1272 template <typename NumberType>
1273 Expression
1274 Expression::substitute(const Expression &symbol,
1275 const NumberType &value) const
1276 {
1277 Assert(SymEngine::is_a<SymEngine::Symbol>(symbol.get_value()),
1278 ExcMessage(
1279 "Substitution with a number that does not represent a symbol."));
1280
1281 types::substitution_map sub_vals;
1282 sub_vals[symbol] = Expression(value);
1283 return substitute(sub_vals);
1284 }
1285
1286
1287 template <typename ReturnType>
1288 ReturnType
1289 Expression::substitute_and_evaluate(
1290 const types::substitution_map &substitution_values) const
1291 {
1292 return static_cast<ReturnType>(substitute(substitution_values));
1293 }
1294
1295
1296 template <typename ReturnType>
1297 ReturnType
1298 Expression::substitute_and_evaluate(
1299 const SymEngine::map_basic_basic &substitution_values) const
1300 {
1301 return static_cast<ReturnType>(substitute(substitution_values));
1302 }
1303
1304
1305 template <typename NumberType>
1306 Expression &
1307 Expression::operator=(const NumberType &rhs)
1308 {
1309 *this = Expression(rhs);
1310 return *this;
1311 }
1312
1313
1314 template <typename NumberType>
1315 Expression &
1316 Expression::operator+=(const NumberType &rhs)
1317 {
1318 *this = Expression(SymEngine::add(get_RCP(), Expression(rhs).get_RCP()));
1319 return *this;
1320 }
1321
1322
1323 template <typename NumberType>
1324 Expression &
1325 Expression::operator-=(const NumberType &rhs)
1326 {
1327 *this = Expression(SymEngine::sub(get_RCP(), Expression(rhs).get_RCP()));
1328 return *this;
1329 }
1330
1331
1332 template <typename NumberType>
1333 Expression &
1334 Expression::operator*=(const NumberType &rhs)
1335 {
1336 *this = Expression(SymEngine::mul(get_RCP(), Expression(rhs).get_RCP()));
1337 return *this;
1338 }
1339
1340
1341 template <typename NumberType>
1342 Expression &
1343 Expression::operator/=(const NumberType &rhs)
1344 {
1345 *this = Expression(SymEngine::div(get_RCP(), Expression(rhs).get_RCP()));
1346 return *this;
1347 }
1348
1349
1350 template <typename ResultType>
1351 Expression::operator ResultType() const
1352 {
1353 return static_cast<ResultType>(get_expression());
1354 }
1355
1356 } // namespace SD
1357} // namespace Differentiation
1358
1359
1360// Specializations from numbers.h
1361// These are required in order to make the SD::Expression class a compatible
1362// number for the Tensor and SymmetricTensor classes
1363
1364// Forward declarations:
1365template <int rank_, int dim, typename Number>
1366class Tensor;
1367
1368namespace internal
1369{
1370 template <>
1371 struct NumberType<Differentiation::SD::Expression>
1372 {
1373 static const Differentiation::SD::Expression &
1375 {
1376 return t;
1377 }
1378
1379 template <
1380 typename T,
1381 typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1383 value(const T &t)
1384 {
1386 }
1387
1388 template <
1389 typename T,
1390 typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1392 value(T &&t)
1393 {
1395 }
1396
1397 template <
1398 typename T,
1399 typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1401 value(const std::complex<T> &t)
1402 {
1404 }
1405
1406 template <typename T, int dim>
1408 value(const Tensor<0, dim, T> &t)
1409 {
1410 return Differentiation::SD::Expression(static_cast<T>(t));
1411 }
1412
1413 template <typename T, int dim>
1415 value(const Tensor<0, dim, std::complex<T>> &t)
1416 {
1417 return Differentiation::SD::Expression(static_cast<std::complex<T>>(t));
1418 }
1419 };
1420} // namespace internal
1421
1422
1423namespace numbers
1424{
1425 template <>
1426 inline bool
1428 {
1429 if (SymEngine::is_a_Number(value.get_value()))
1430 {
1431 const SymEngine::RCP<const SymEngine::Number> number_rcp =
1432 SymEngine::rcp_static_cast<const SymEngine::Number>(value.get_RCP());
1433 return number_rcp->is_zero();
1434 }
1435
1436 return false;
1437 }
1438
1439 template <>
1440 inline bool
1442 const Differentiation::SD::Expression &value_2)
1443 {
1444 return (value_1.get_value().__cmp__(value_2.get_value()) == 0);
1445 }
1446
1447 template <>
1448 inline bool
1450 const Differentiation::SD::Expression &value_2)
1451 {
1452 return (value_1.get_value().__cmp__(value_2.get_value()) == -1);
1453 }
1454} // namespace numbers
1455
1456
1457# endif // DOXYGEN
1458
1460
1461#endif // DEAL_II_WITH_SYMENGINE
1462
1463#endif // dealii_differentiation_sd_symengine_number_types_h
Definition: tensor.h:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Expression & operator/=(const Expression &rhs)
Expression & parse(const std::string &expression)
Expression & operator-=(const NumberType &rhs)
Expression operator*(Expression lhs, const Expression &rhs)
Expression operator<(const Expression &lhs, const Expression &rhs)
Expression operator-(Expression lhs, const Expression &rhs)
Expression substitute(const Expression &symbol, const NumberType &value) const
Expression & operator=(const Expression &rhs)
Expression operator+(Expression lhs, const Expression &rhs)
Expression(const NumberType &value)
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
ReturnType substitute_and_evaluate(const types::substitution_map &substitution_values) const
Expression & operator/=(const NumberType &rhs)
Expression(const NumberType &numerator, const NumberType &denominator)
Expression(Expression &&rhs)=default
Expression & operator*=(const Expression &rhs)
Expression operator!(const Expression &expression)
Expression operator||(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcSymEngineParserError(std::string arg1)
Expression operator^(const Expression &lhs, const Expression &rhs)
void serialize(Archive &archive, const unsigned int version)
#define Assert(cond, exc)
Definition: exceptions.h:1473
Expression & operator+=(const NumberType &rhs)
Expression operator>=(const Expression &lhs, const Expression &rhs)
std::ostream & print(std::ostream &stream) const
Expression substitute(const types::substitution_map &substitution_values) const
Expression operator!=(const Expression &lhs, const Expression &rhs)
Expression & operator=(const NumberType &rhs)
Expression operator|(const Expression &lhs, const Expression &rhs)
void save(Archive &archive, const unsigned int version) const
Expression(const Expression &rhs)=default
void save(std::ostream &stream) const
Expression & operator*=(const NumberType &rhs)
ReturnType substitute_and_evaluate(const SymEngine::map_basic_basic &substitution_values) const
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 &rhs)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
Expression operator&&(const Expression &lhs, const Expression &rhs)
const SymEngine::Basic & get_value() const
const SymEngine::Expression & get_expression() const
std::ostream & operator<<(std::ostream &stream, const Expression &expression)
Expression operator==(const Expression &lhs, const Expression &rhs)
Expression & operator+=(const Expression &rhs)
void load(Archive &archive, const unsigned int version)
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression(const std::complex< NumberType > &value)
Expression differentiate(const Expression &symbol) const
std::istream & operator>>(std::istream &stream, Expression &expression)
std::vector< SD::Expression > symbol_vector
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:943
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:927
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:951