Reference documentation for deal.II version Git a7d922cd5a 2021-01-15 21:19:32 -0700
\(\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.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 - 2020 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
48 # include <deal.II/base/exceptions.h>
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 
67 namespace 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
483 
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 
619  Expression
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 
673  Expression
674  differentiate(const Expression &symbol) const;
675 
680  Expression
682  const SymEngine::RCP<const SymEngine::Symbol> &symbol) const;
683 
688  Expression
689  differentiate(const SymEngine::RCP<const SymEngine::Basic> &symbol) const;
690 
692 
697 
710  Expression
711  substitute(const types::substitution_map &substitution_values) const;
712 
729  Expression
730  substitute(const SymEngine::map_basic_basic &substitution_values) const;
731 
741  Expression
742  substitute(const Expression &symbol, const Expression &value) const;
743 
753  template <typename NumberType>
754  Expression
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 &
864  get_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 
933  Expression
934  operator==(const Expression &lhs, const Expression &rhs);
935 
941  Expression
942  operator!=(const Expression &lhs, const Expression &rhs);
943 
949  Expression
950  operator<(const Expression &lhs, const Expression &rhs);
951 
957  Expression
958  operator>(const Expression &lhs, const Expression &rhs);
959 
965  Expression
966  operator<=(const Expression &lhs, const Expression &rhs);
967 
973  Expression
974  operator>=(const Expression &lhs, const Expression &rhs);
975 
977 
982 
989  Expression operator!(const Expression &expression);
990 
997  Expression operator&(const Expression &lhs, const Expression &rhs);
998 
1005  Expression
1006  operator|(const Expression &lhs, const Expression &rhs);
1007 
1014  Expression
1015  operator^(const Expression &lhs, const Expression &rhs);
1016 
1024  Expression
1025  operator&&(const Expression &lhs, const Expression &rhs);
1026 
1034  Expression
1035  operator||(const Expression &lhs, const Expression &rhs);
1036 
1038 
1043 
1049  Expression
1050  operator+(Expression lhs, const Expression &rhs);
1051 
1057  Expression
1058  operator-(Expression lhs, const Expression &rhs);
1059 
1065  Expression operator*(Expression lhs, const Expression &rhs);
1066 
1072  Expression
1073  operator/(Expression lhs, const Expression &rhs);
1074 
1083  template <typename NumberType,
1084  typename = typename std::enable_if<
1085  std::is_constructible<Expression, NumberType>::value>::type>
1086  inline Expression
1087  operator+(const NumberType &lhs, const Expression &rhs)
1088  {
1089  return Expression(lhs) + rhs;
1090  }
1091 
1100  template <typename NumberType,
1101  typename = typename std::enable_if<
1102  std::is_constructible<Expression, NumberType>::value>::type>
1103  inline Expression
1104  operator+(const Expression &lhs, const NumberType &rhs)
1105  {
1106  return lhs + Expression(rhs);
1107  }
1108 
1117  template <typename NumberType,
1118  typename = typename std::enable_if<
1119  std::is_constructible<Expression, NumberType>::value>::type>
1120  inline Expression
1121  operator-(const NumberType &lhs, const Expression &rhs)
1122  {
1123  return Expression(lhs) - rhs;
1124  }
1125 
1134  template <typename NumberType,
1135  typename = typename std::enable_if<
1136  std::is_constructible<Expression, NumberType>::value>::type>
1137  inline Expression
1138  operator-(const Expression &lhs, const NumberType &rhs)
1139  {
1140  return lhs - Expression(rhs);
1141  }
1142 
1151  template <typename NumberType,
1152  typename = typename std::enable_if<
1153  std::is_constructible<Expression, NumberType>::value>::type>
1154  inline Expression operator*(const NumberType &lhs, const Expression &rhs)
1155  {
1156  return Expression(lhs) * rhs;
1157  }
1158 
1167  template <typename NumberType,
1168  typename = typename std::enable_if<
1169  std::is_constructible<Expression, NumberType>::value>::type>
1170  inline Expression operator*(const Expression &lhs, const NumberType &rhs)
1171  {
1172  return lhs * Expression(rhs);
1173  }
1174 
1183  template <typename NumberType,
1184  typename = typename std::enable_if<
1185  std::is_constructible<Expression, NumberType>::value>::type>
1186  inline Expression
1187  operator/(const NumberType &lhs, const Expression &rhs)
1188  {
1189  return Expression(lhs) / rhs;
1190  }
1191 
1200  template <typename NumberType,
1201  typename = typename std::enable_if<
1202  std::is_constructible<Expression, NumberType>::value>::type>
1203  inline Expression
1204  operator/(const Expression &lhs, const NumberType &rhs)
1205  {
1206  return lhs / Expression(rhs);
1207  }
1208 
1210 
1211  } // namespace SD
1212 } // namespace Differentiation
1213 
1214 
1215 /* ---------------------- inline and template functions -------------------- */
1216 
1217 
1218 # ifndef DOXYGEN
1219 
1220 
1221 namespace Differentiation
1222 {
1223  namespace SD
1224  {
1225  template <typename NumberType, typename>
1226  Expression::Expression(const NumberType &value)
1227  : expression(value)
1228  {}
1229 
1230 
1231  template <typename NumberType, typename>
1232  Expression::Expression(const std::complex<NumberType> &value)
1233  : expression(value)
1234  {}
1235 
1236 
1237  template <typename NumberType, typename>
1238  Expression::Expression(const NumberType &numerator,
1239  const NumberType &denominator)
1240  : expression(
1241  SymEngine::Rational::from_two_ints(*SymEngine::integer(numerator),
1242  *SymEngine::integer(denominator)))
1243  {}
1244 
1245 
1246  template <class Archive>
1247  void
1248  Expression::save(Archive &ar, const unsigned int /*version*/) const
1249  {
1250  std::stringstream sstream;
1251  sstream << *this;
1252  const std::string expr = sstream.str();
1253  ar & expr;
1254  }
1255 
1256 
1257  template <class Archive>
1258  void
1259  Expression::load(Archive &ar, const unsigned int /*version*/)
1260  {
1261  std::string expr;
1262  ar & expr;
1263  parse(expr);
1264  }
1265 
1266 
1267  template <typename NumberType>
1268  Expression
1269  Expression::substitute(const Expression &symbol,
1270  const NumberType &value) const
1271  {
1272  Assert(SymEngine::is_a<SymEngine::Symbol>(symbol.get_value()),
1273  ExcMessage(
1274  "Substitution with a number that does not represent a symbol."));
1275 
1276  types::substitution_map sub_vals;
1277  sub_vals[symbol] = Expression(value);
1278  return substitute(sub_vals);
1279  }
1280 
1281 
1282  template <typename ReturnType>
1283  ReturnType
1285  const types::substitution_map &substitution_values) const
1286  {
1287  return static_cast<ReturnType>(substitute(substitution_values));
1288  }
1289 
1290 
1291  template <typename ReturnType>
1292  ReturnType
1294  const SymEngine::map_basic_basic &substitution_values) const
1295  {
1296  return static_cast<ReturnType>(substitute(substitution_values));
1297  }
1298 
1299 
1300  template <typename NumberType>
1301  Expression &
1302  Expression::operator=(const NumberType &rhs)
1303  {
1304  *this = Expression(rhs);
1305  return *this;
1306  }
1307 
1308 
1309  template <typename NumberType>
1310  Expression &
1311  Expression::operator+=(const NumberType &rhs)
1312  {
1313  *this = Expression(SymEngine::add(get_RCP(), Expression(rhs).get_RCP()));
1314  return *this;
1315  }
1316 
1317 
1318  template <typename NumberType>
1319  Expression &
1320  Expression::operator-=(const NumberType &rhs)
1321  {
1322  *this = Expression(SymEngine::sub(get_RCP(), Expression(rhs).get_RCP()));
1323  return *this;
1324  }
1325 
1326 
1327  template <typename NumberType>
1328  Expression &
1329  Expression::operator*=(const NumberType &rhs)
1330  {
1331  *this = Expression(SymEngine::mul(get_RCP(), Expression(rhs).get_RCP()));
1332  return *this;
1333  }
1334 
1335 
1336  template <typename NumberType>
1337  Expression &
1338  Expression::operator/=(const NumberType &rhs)
1339  {
1340  *this = Expression(SymEngine::div(get_RCP(), Expression(rhs).get_RCP()));
1341  return *this;
1342  }
1343 
1344 
1345  template <typename ResultType>
1346  Expression::operator ResultType() const
1347  {
1348  return static_cast<ResultType>(get_expression());
1349  }
1350 
1351  } // namespace SD
1352 } // namespace Differentiation
1353 
1354 
1355 // Specializations from numbers.h
1356 // These are required in order to make the SD::Expression class a compatible
1357 // number for the Tensor and SymmetricTensor classes
1358 
1359 // Forward declarations:
1360 template <int rank_, int dim, typename Number>
1361 class Tensor;
1362 
1363 namespace internal
1364 {
1365  template <>
1366  struct NumberType<Differentiation::SD::Expression>
1367  {
1368  static const Differentiation::SD::Expression &
1369  value(const Differentiation::SD::Expression &t)
1370  {
1371  return t;
1372  }
1373 
1374  template <
1375  typename T,
1376  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1378  value(const T &t)
1379  {
1381  }
1382 
1383  template <
1384  typename T,
1385  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1387  value(T &&t)
1388  {
1390  }
1391 
1392  template <
1393  typename T,
1394  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1396  value(const std::complex<T> &t)
1397  {
1399  }
1400 
1401  template <typename T, int dim>
1403  value(const Tensor<0, dim, T> &t)
1404  {
1405  return Differentiation::SD::Expression(static_cast<T>(t));
1406  }
1407 
1408  template <typename T, int dim>
1410  value(const Tensor<0, dim, std::complex<T>> &t)
1411  {
1412  return Differentiation::SD::Expression(static_cast<std::complex<T>>(t));
1413  }
1414  };
1415 } // namespace internal
1416 
1417 
1418 namespace numbers
1419 {
1420  template <>
1421  inline bool
1423  {
1424  if (SymEngine::is_a_Number(value.get_value()))
1425  {
1426  const SymEngine::RCP<const SymEngine::Number> number_rcp =
1427  SymEngine::rcp_static_cast<const SymEngine::Number>(value.get_RCP());
1428  return number_rcp->is_zero();
1429  }
1430 
1431  return false;
1432  }
1433 
1434  template <>
1435  inline bool
1437  const Differentiation::SD::Expression &value_2)
1438  {
1439  return (value_1.get_value().__cmp__(value_2.get_value()) == 0);
1440  }
1441 
1442  template <>
1443  inline bool
1445  const Differentiation::SD::Expression &value_2)
1446  {
1447  return (value_1.get_value().__cmp__(value_2.get_value()) == -1);
1448  }
1449 } // namespace numbers
1450 
1451 
1452 # endif // DOXYGEN
1453 
1455 
1456 #endif // DEAL_II_WITH_SYMENGINE
1457 
1458 #endif // dealii_differentiation_sd_symengine_number_types_h
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)
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:915
Expression operator>=(const Expression &lhs, const Expression &rhs)
Expression & parse(const std::string &expression)
ReturnType substitute_and_evaluate(const types::substitution_map &substitution_values) const
void serialize(Archive &archive, const unsigned int version)
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)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
static const char T
Expression operator<=(const Expression &lhs, const Expression &rhs)
#define Assert(cond, exc)
Definition: exceptions.h:1466
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
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:380
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
std::vector< SD::Expression > symbol_vector
Expression operator*(Expression lhs, const Expression &rhs)
Expression operator!(const Expression &expression)
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:939
Definition: tensor.h:449
Expression operator==(const Expression &lhs, const Expression &rhs)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:379
Expression operator||(const Expression &lhs, const Expression &rhs)
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:931
Expression operator^(const Expression &lhs, const Expression &rhs)
Expression & operator-=(const Expression &rhs)
Expression & operator+=(const Expression &rhs)
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
Expression operator!=(const Expression &lhs, const Expression &rhs)
Expression operator &(const Expression &lhs, const Expression &rhs)