Reference documentation for deal.II version Git 662679e6de 2020-10-27 13:02:35 -0600
\(\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 
450  void
451  load(std::istream &stream);
452 
460  template <class Archive>
461  void
462  save(Archive &archive, const unsigned int version) const;
463 
476  template <class Archive>
477  void
478  load(Archive &archive, const unsigned int version);
479 
480 # ifdef DOXYGEN
481 
494  template <class Archive>
495  void
496  serialize(Archive &archive, const unsigned int version);
497 # else
498  // This macro defines the serialize() method that is compatible with
499  // the templated save() and load() method that have been implemented.
500  BOOST_SERIALIZATION_SPLIT_MEMBER()
501 # endif
502 
504 
509 
513  const SymEngine::Expression &
514  get_expression() const;
515 
520  const SymEngine::Basic &
521  get_value() const;
522 
527  const SymEngine::RCP<const SymEngine::Basic> &
528  get_RCP() const;
529 
531 
536 
543  Expression &
544  operator=(const Expression &rhs);
545 
552  Expression &
553  operator=(Expression &&rhs) noexcept;
554 
561  Expression &
562  operator+=(const Expression &rhs);
563 
570  Expression &
571  operator-=(const Expression &rhs);
572 
579  Expression &
580  operator*=(const Expression &rhs);
581 
588  Expression &
589  operator/=(const Expression &rhs);
590 
592 
597 
604  template <typename NumberType>
605  Expression &
606  operator=(const NumberType &rhs);
607 
616  Expression
617  operator-() const;
618 
625  template <typename NumberType>
626  Expression &
627  operator+=(const NumberType &rhs);
628 
635  template <typename NumberType>
636  Expression &
637  operator-=(const NumberType &rhs);
638 
645  template <typename NumberType>
646  Expression &
647  operator*=(const NumberType &rhs);
648 
655  template <typename NumberType>
656  Expression &
657  operator/=(const NumberType &rhs);
658 
660 
665 
670  Expression
671  differentiate(const Expression &symbol) const;
672 
677  Expression
679  const SymEngine::RCP<const SymEngine::Symbol> &symbol) const;
680 
685  Expression
686  differentiate(const SymEngine::RCP<const SymEngine::Basic> &symbol) const;
687 
689 
694 
707  Expression
708  substitute(const types::substitution_map &substitution_values) const;
709 
726  Expression
727  substitute(const SymEngine::map_basic_basic &substitution_values) const;
728 
738  Expression
739  substitute(const Expression &symbol, const Expression &value) const;
740 
750  template <typename NumberType>
751  Expression
752  substitute(const Expression &symbol, const NumberType &value) const;
753 
761  template <typename ReturnType>
762  ReturnType
764  const types::substitution_map &substitution_values) const;
765 
777  template <typename ReturnType>
778  ReturnType
780  const SymEngine::map_basic_basic &substitution_values) const;
781 
783 
788 
839  template <typename ResultType>
840  explicit operator ResultType() const;
841 
846  explicit operator const SymEngine::Expression &() const;
847 
852  operator const SymEngine::RCP<const SymEngine::Basic> &() const;
853 
855 
856  protected:
860  SymEngine::Expression &
861  get_expression();
862 
863  private:
868  SymEngine::Expression expression;
869  };
870 
875 
881  template <>
882  struct is_symengine_number<Expression> : std::true_type
883  {};
884 
885 
891  template <>
892  struct is_sd_number<Expression> : std::true_type
893  {};
894 
896 
901 
907  std::ostream &
908  operator<<(std::ostream &stream, const Expression &expression);
909 
915  std::istream &
916  operator>>(std::istream &stream, Expression &expression);
917 
919 
924 
930  Expression
931  operator==(const Expression &lhs, const Expression &rhs);
932 
938  Expression
939  operator!=(const Expression &lhs, const Expression &rhs);
940 
946  Expression
947  operator<(const Expression &lhs, const Expression &rhs);
948 
954  Expression
955  operator>(const Expression &lhs, const Expression &rhs);
956 
962  Expression
963  operator<=(const Expression &lhs, const Expression &rhs);
964 
970  Expression
971  operator>=(const Expression &lhs, const Expression &rhs);
972 
974 
979 
986  Expression operator!(const Expression &expression);
987 
994  Expression operator&(const Expression &lhs, const Expression &rhs);
995 
1002  Expression
1003  operator|(const Expression &lhs, const Expression &rhs);
1004 
1011  Expression
1012  operator^(const Expression &lhs, const Expression &rhs);
1013 
1021  Expression
1022  operator&&(const Expression &lhs, const Expression &rhs);
1023 
1031  Expression
1032  operator||(const Expression &lhs, const Expression &rhs);
1033 
1035 
1040 
1046  Expression
1047  operator+(Expression lhs, const Expression &rhs);
1048 
1054  Expression
1055  operator-(Expression lhs, const Expression &rhs);
1056 
1062  Expression operator*(Expression lhs, const Expression &rhs);
1063 
1069  Expression
1070  operator/(Expression lhs, const Expression &rhs);
1071 
1080  template <typename NumberType,
1081  typename = typename std::enable_if<
1082  std::is_constructible<Expression, NumberType>::value>::type>
1083  inline Expression
1084  operator+(const NumberType &lhs, const Expression &rhs)
1085  {
1086  return Expression(lhs) + rhs;
1087  }
1088 
1097  template <typename NumberType,
1098  typename = typename std::enable_if<
1099  std::is_constructible<Expression, NumberType>::value>::type>
1100  inline Expression
1101  operator+(const Expression &lhs, const NumberType &rhs)
1102  {
1103  return lhs + Expression(rhs);
1104  }
1105 
1114  template <typename NumberType,
1115  typename = typename std::enable_if<
1116  std::is_constructible<Expression, NumberType>::value>::type>
1117  inline Expression
1118  operator-(const NumberType &lhs, const Expression &rhs)
1119  {
1120  return Expression(lhs) - rhs;
1121  }
1122 
1131  template <typename NumberType,
1132  typename = typename std::enable_if<
1133  std::is_constructible<Expression, NumberType>::value>::type>
1134  inline Expression
1135  operator-(const Expression &lhs, const NumberType &rhs)
1136  {
1137  return lhs - Expression(rhs);
1138  }
1139 
1148  template <typename NumberType,
1149  typename = typename std::enable_if<
1150  std::is_constructible<Expression, NumberType>::value>::type>
1151  inline Expression operator*(const NumberType &lhs, const Expression &rhs)
1152  {
1153  return Expression(lhs) * rhs;
1154  }
1155 
1164  template <typename NumberType,
1165  typename = typename std::enable_if<
1166  std::is_constructible<Expression, NumberType>::value>::type>
1167  inline Expression operator*(const Expression &lhs, const NumberType &rhs)
1168  {
1169  return lhs * Expression(rhs);
1170  }
1171 
1180  template <typename NumberType,
1181  typename = typename std::enable_if<
1182  std::is_constructible<Expression, NumberType>::value>::type>
1183  inline Expression
1184  operator/(const NumberType &lhs, const Expression &rhs)
1185  {
1186  return Expression(lhs) / rhs;
1187  }
1188 
1197  template <typename NumberType,
1198  typename = typename std::enable_if<
1199  std::is_constructible<Expression, NumberType>::value>::type>
1200  inline Expression
1201  operator/(const Expression &lhs, const NumberType &rhs)
1202  {
1203  return lhs / Expression(rhs);
1204  }
1205 
1207 
1208  } // namespace SD
1209 } // namespace Differentiation
1210 
1211 
1212 /* ---------------------- inline and template functions -------------------- */
1213 
1214 
1215 # ifndef DOXYGEN
1216 
1217 
1218 namespace Differentiation
1219 {
1220  namespace SD
1221  {
1222  template <typename NumberType, typename>
1223  Expression::Expression(const NumberType &value)
1224  : expression(value)
1225  {}
1226 
1227 
1228  template <typename NumberType, typename>
1229  Expression::Expression(const std::complex<NumberType> &value)
1230  : expression(value)
1231  {}
1232 
1233 
1234  template <typename NumberType, typename>
1235  Expression::Expression(const NumberType &numerator,
1236  const NumberType &denominator)
1237  : expression(
1238  SymEngine::Rational::from_two_ints(*SymEngine::integer(numerator),
1239  *SymEngine::integer(denominator)))
1240  {}
1241 
1242 
1243  template <class Archive>
1244  void
1245  Expression::save(Archive &ar, const unsigned int /*version*/) const
1246  {
1247  std::stringstream sstream;
1248  sstream << *this;
1249  const std::string expr = sstream.str();
1250  ar & expr;
1251  }
1252 
1253 
1254  template <class Archive>
1255  void
1256  Expression::load(Archive &ar, const unsigned int /*version*/)
1257  {
1258  std::string expr;
1259  ar & expr;
1260  parse(expr);
1261  }
1262 
1263 
1264  template <typename NumberType>
1265  Expression
1266  Expression::substitute(const Expression &symbol,
1267  const NumberType &value) const
1268  {
1269  Assert(SymEngine::is_a<SymEngine::Symbol>(symbol.get_value()),
1270  ExcMessage(
1271  "Substitution with a number that does not represent a symbol."));
1272 
1273  types::substitution_map sub_vals;
1274  sub_vals[symbol] = Expression(value);
1275  return substitute(sub_vals);
1276  }
1277 
1278 
1279  template <typename ReturnType>
1280  ReturnType
1282  const types::substitution_map &substitution_values) const
1283  {
1284  return static_cast<ReturnType>(substitute(substitution_values));
1285  }
1286 
1287 
1288  template <typename ReturnType>
1289  ReturnType
1291  const SymEngine::map_basic_basic &substitution_values) const
1292  {
1293  return static_cast<ReturnType>(substitute(substitution_values));
1294  }
1295 
1296 
1297  template <typename NumberType>
1298  Expression &
1299  Expression::operator=(const NumberType &rhs)
1300  {
1301  *this = Expression(rhs);
1302  return *this;
1303  }
1304 
1305 
1306  template <typename NumberType>
1307  Expression &
1308  Expression::operator+=(const NumberType &rhs)
1309  {
1310  *this = Expression(SymEngine::add(get_RCP(), Expression(rhs).get_RCP()));
1311  return *this;
1312  }
1313 
1314 
1315  template <typename NumberType>
1316  Expression &
1317  Expression::operator-=(const NumberType &rhs)
1318  {
1319  *this = Expression(SymEngine::sub(get_RCP(), Expression(rhs).get_RCP()));
1320  return *this;
1321  }
1322 
1323 
1324  template <typename NumberType>
1325  Expression &
1326  Expression::operator*=(const NumberType &rhs)
1327  {
1328  *this = Expression(SymEngine::mul(get_RCP(), Expression(rhs).get_RCP()));
1329  return *this;
1330  }
1331 
1332 
1333  template <typename NumberType>
1334  Expression &
1335  Expression::operator/=(const NumberType &rhs)
1336  {
1337  *this = Expression(SymEngine::div(get_RCP(), Expression(rhs).get_RCP()));
1338  return *this;
1339  }
1340 
1341 
1342  template <typename ResultType>
1343  Expression::operator ResultType() const
1344  {
1345  return static_cast<ResultType>(get_expression());
1346  }
1347 
1348  } // namespace SD
1349 } // namespace Differentiation
1350 
1351 
1352 // Specializations from numbers.h
1353 // These are required in order to make the SD::Expression class a compatible
1354 // number for the Tensor and SymmetricTensor classes
1355 
1356 // Forward declarations:
1357 template <int rank_, int dim, typename Number>
1358 class Tensor;
1359 
1360 namespace internal
1361 {
1362  template <>
1363  struct NumberType<Differentiation::SD::Expression>
1364  {
1365  static const Differentiation::SD::Expression &
1366  value(const Differentiation::SD::Expression &t)
1367  {
1368  return t;
1369  }
1370 
1371  template <
1372  typename T,
1373  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1375  value(const T &t)
1376  {
1378  }
1379 
1380  template <
1381  typename T,
1382  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1384  value(T &&t)
1385  {
1387  }
1388 
1389  template <
1390  typename T,
1391  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1393  value(const std::complex<T> &t)
1394  {
1396  }
1397 
1398  template <typename T, int dim>
1400  value(const Tensor<0, dim, T> &t)
1401  {
1402  return Differentiation::SD::Expression(static_cast<T>(t));
1403  }
1404 
1405  template <typename T, int dim>
1407  value(const Tensor<0, dim, std::complex<T>> &t)
1408  {
1409  return Differentiation::SD::Expression(static_cast<std::complex<T>>(t));
1410  }
1411  };
1412 } // namespace internal
1413 
1414 
1415 namespace numbers
1416 {
1417  template <>
1418  inline bool
1420  {
1421  if (SymEngine::is_a_Number(value.get_value()))
1422  {
1423  const SymEngine::RCP<const SymEngine::Number> number_rcp =
1424  SymEngine::rcp_static_cast<const SymEngine::Number>(value.get_RCP());
1425  return number_rcp->is_zero();
1426  }
1427 
1428  return false;
1429  }
1430 
1431  template <>
1432  inline bool
1434  const Differentiation::SD::Expression &value_2)
1435  {
1436  return (value_1.get_value().__cmp__(value_2.get_value()) == 0);
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()) == -1);
1445  }
1446 } // namespace numbers
1447 
1448 
1449 # endif // DOXYGEN
1450 
1452 
1453 #endif // DEAL_II_WITH_SYMENGINE
1454 
1455 #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:1423
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:369
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:448
Expression operator==(const Expression &lhs, const Expression &rhs)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
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)