Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
symengine_number_types.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 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 
51 # include <deal.II/differentiation/sd/symengine_number_traits.h>
52 # include <deal.II/differentiation/sd/symengine_types.h>
53 
54 # include <symengine/derivative.h>
55 
56 # include <algorithm>
57 # include <memory>
58 # include <sstream>
59 # include <type_traits>
60 # include <utility>
61 # include <vector>
62 
63 DEAL_II_NAMESPACE_OPEN
64 
65 namespace Differentiation
66 {
67  namespace SD
68  {
80  std::string,
81  << "The string '" << arg1
82  << "' could not be parsed successfully. Are you sure that (1) it "
83  << "consists of legitimate operations and syntax, and (2) you've "
84  << "previously declared all symbolic variables that are present "
85  << "in the expression?");
86 
88 
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 
466  template <class Archive>
467  void
468  serialize(Archive &archive, const unsigned int version);
469 
471 
476 
480  const SymEngine::Expression &
481  get_expression() const;
482 
487  const SymEngine::Basic &
488  get_value() const;
489 
494  const SymEngine::RCP<const SymEngine::Basic> &
495  get_RCP() const;
496 
498 
503 
510  Expression &
511  operator=(const Expression &rhs);
512 
519  Expression &
520  operator=(Expression &&rhs);
521 
528  Expression &
529  operator+=(const Expression &rhs);
530 
537  Expression &
538  operator-=(const Expression &rhs);
539 
546  Expression &
547  operator*=(const Expression &rhs);
548 
555  Expression &
556  operator/=(const Expression &rhs);
557 
559 
564 
571  template <typename NumberType>
572  Expression &
573  operator=(const NumberType &rhs);
574 
583  Expression
584  operator-() const;
585 
592  template <typename NumberType>
593  Expression &
594  operator+=(const NumberType &rhs);
595 
602  template <typename NumberType>
603  Expression &
604  operator-=(const NumberType &rhs);
605 
612  template <typename NumberType>
613  Expression &
614  operator*=(const NumberType &rhs);
615 
622  template <typename NumberType>
623  Expression &
624  operator/=(const NumberType &rhs);
625 
627 
632 
637  Expression
638  differentiate(const Expression &symbol) const;
639 
644  Expression
646  const SymEngine::RCP<const SymEngine::Symbol> &symbol) const;
647 
652  Expression
653  differentiate(const SymEngine::RCP<const SymEngine::Basic> &symbol) const;
654 
656 
661 
674  Expression
675  substitute(const types::substitution_map &substitution_values) const;
676 
686  Expression
687  substitute(const Expression &symbol, const Expression &value) const;
688 
698  template <typename NumberType>
699  Expression
700  substitute(const Expression &symbol, const NumberType &value) const;
701 
709  template <typename ReturnType>
710  ReturnType
712  const types::substitution_map &substitution_values) const;
713 
715 
720 
771  template <typename ResultType>
772  explicit operator ResultType() const;
773 
778  explicit operator const SymEngine::Expression &() const;
779 
784  operator const SymEngine::RCP<const SymEngine::Basic> &() const;
785 
787 
788  protected:
792  SymEngine::Expression &
793  get_expression();
794 
795  private:
800  SymEngine::Expression expression;
801  };
802 
807 
815  template <>
816  struct is_symengine_number<Expression> : std::true_type
817  {};
818 
819 
827  template <>
828  struct is_sd_number<Expression> : std::true_type
829  {};
830 
832 
837 
843  std::ostream &
844  operator<<(std::ostream &stream, const Expression &expression);
845 
851  std::istream &
852  operator>>(std::istream &stream, Expression &expression);
853 
855 
860 
866  Expression
867  operator==(const Expression &lhs, const Expression &rhs);
868 
874  Expression
875  operator!=(const Expression &lhs, const Expression &rhs);
876 
882  Expression
883  operator<(const Expression &lhs, const Expression &rhs);
884 
890  Expression
891  operator>(const Expression &lhs, const Expression &rhs);
892 
898  Expression
899  operator<=(const Expression &lhs, const Expression &rhs);
900 
906  Expression
907  operator>=(const Expression &lhs, const Expression &rhs);
908 
910 
915 
922  Expression operator!(const Expression &expression);
923 
930  Expression operator&(const Expression &lhs, const Expression &rhs);
931 
938  Expression
939  operator|(const Expression &lhs, const Expression &rhs);
940 
947  Expression
948  operator^(const Expression &lhs, const Expression &rhs);
949 
957  Expression
958  operator&&(const Expression &lhs, const Expression &rhs);
959 
967  Expression
968  operator||(const Expression &lhs, const Expression &rhs);
969 
971 
976 
982  Expression
983  operator+(Expression lhs, const Expression &rhs);
984 
990  Expression
991  operator-(Expression lhs, const Expression &rhs);
992 
998  Expression operator*(Expression lhs, const Expression &rhs);
999 
1005  Expression
1006  operator/(Expression lhs, const Expression &rhs);
1007 
1016  template <typename NumberType,
1017  typename = typename std::enable_if<
1018  std::is_constructible<Expression, NumberType>::value>::type>
1019  inline Expression
1020  operator+(const NumberType &lhs, const Expression &rhs)
1021  {
1022  return Expression(lhs) + rhs;
1023  }
1024 
1033  template <typename NumberType,
1034  typename = typename std::enable_if<
1035  std::is_constructible<Expression, NumberType>::value>::type>
1036  inline Expression
1037  operator+(const Expression &lhs, const NumberType &rhs)
1038  {
1039  return lhs + Expression(rhs);
1040  }
1041 
1050  template <typename NumberType,
1051  typename = typename std::enable_if<
1052  std::is_constructible<Expression, NumberType>::value>::type>
1053  inline Expression
1054  operator-(const NumberType &lhs, const Expression &rhs)
1055  {
1056  return Expression(lhs) - rhs;
1057  }
1058 
1067  template <typename NumberType,
1068  typename = typename std::enable_if<
1069  std::is_constructible<Expression, NumberType>::value>::type>
1070  inline Expression
1071  operator-(const Expression &lhs, const NumberType &rhs)
1072  {
1073  return lhs - Expression(rhs);
1074  }
1075 
1084  template <typename NumberType,
1085  typename = typename std::enable_if<
1086  std::is_constructible<Expression, NumberType>::value>::type>
1087  inline Expression 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 operator*(const Expression &lhs, const NumberType &rhs)
1104  {
1105  return lhs * Expression(rhs);
1106  }
1107 
1116  template <typename NumberType,
1117  typename = typename std::enable_if<
1118  std::is_constructible<Expression, NumberType>::value>::type>
1119  inline Expression
1120  operator/(const NumberType &lhs, const Expression &rhs)
1121  {
1122  return Expression(lhs) / rhs;
1123  }
1124 
1133  template <typename NumberType,
1134  typename = typename std::enable_if<
1135  std::is_constructible<Expression, NumberType>::value>::type>
1136  inline Expression
1137  operator/(const Expression &lhs, const NumberType &rhs)
1138  {
1139  return lhs / Expression(rhs);
1140  }
1141 
1143 
1144  } // namespace SD
1145 } // namespace Differentiation
1146 
1147 
1148 /* ---------------------- inline and template functions -------------------- */
1149 
1150 
1151 # ifndef DOXYGEN
1152 
1153 
1154 namespace Differentiation
1155 {
1156  namespace SD
1157  {
1158  template <typename NumberType, typename>
1159  Expression::Expression(const NumberType &value)
1160  : expression(value)
1161  {}
1162 
1163 
1164  template <typename NumberType, typename>
1165  Expression::Expression(const std::complex<NumberType> &value)
1166  : expression(value)
1167  {}
1168 
1169 
1170  template <typename NumberType, typename>
1171  Expression::Expression(const NumberType &numerator,
1172  const NumberType &denominator)
1173  : expression(
1174  SymEngine::Rational::from_two_ints(*SymEngine::integer(numerator),
1175  *SymEngine::integer(denominator)))
1176  {}
1177 
1178 
1179  template <class Archive>
1180  void
1181  Expression::serialize(Archive &ar, const unsigned int /*version*/)
1182  {
1183  // This is a bit tricky... The SymEngine expression class is not
1184  // serializable, and we're not sure if we're writing out or
1185  // reading in here. So what we'll do is try to synchronise the
1186  // current data with the stream, predicting if we're likely trying
1187  // to write out data. We do this by checking if we're currently working
1188  // with an object that has a trivial state. If not, then we're likely
1189  // to output this object's contents.
1190  const Expression default_constructed;
1191  const bool likely_writing_out =
1192  (get_RCP()->__eq__(*default_constructed.get_RCP()) == false);
1193 
1194  std::stringstream sstream;
1195  std::string expr = default_constructed.get_RCP()->__str__();
1196 
1197  // Output
1198  if (likely_writing_out)
1199  {
1200  sstream << *this;
1201  expr = sstream.str();
1202  }
1203 
1204  // Serialise
1205  ar &expr;
1206 
1207  // Input
1208  if (!likely_writing_out)
1209  {
1210  sstream.clear();
1211  sstream << expr;
1212  sstream >> *this;
1213 
1214  parse(expr);
1215  }
1216  }
1217 
1218 
1219  template <typename NumberType>
1220  Expression
1221  Expression::substitute(const Expression &symbol,
1222  const NumberType &value) const
1223  {
1224  Assert(SymEngine::is_a<SymEngine::Symbol>(symbol.get_value()),
1225  ExcMessage(
1226  "Substitution with a number that does not represent a symbol."));
1227 
1228  types::substitution_map sub_vals;
1229  sub_vals[symbol] = Expression(value);
1230  return substitute(sub_vals);
1231  }
1232 
1233 
1234  template <typename ReturnType>
1235  ReturnType
1236  Expression::substitute_and_evaluate(
1237  const types::substitution_map &substitution_values) const
1238  {
1239  return static_cast<ReturnType>(substitute(substitution_values));
1240  }
1241 
1242 
1243  template <typename NumberType>
1244  Expression &
1245  Expression::operator=(const NumberType &rhs)
1246  {
1247  *this = Expression(rhs);
1248  return *this;
1249  }
1250 
1251 
1252  template <typename NumberType>
1253  Expression &
1254  Expression::operator+=(const NumberType &rhs)
1255  {
1256  *this = Expression(SymEngine::add(get_RCP(), Expression(rhs).get_RCP()));
1257  return *this;
1258  }
1259 
1260 
1261  template <typename NumberType>
1262  Expression &
1263  Expression::operator-=(const NumberType &rhs)
1264  {
1265  *this = Expression(SymEngine::sub(get_RCP(), Expression(rhs).get_RCP()));
1266  return *this;
1267  }
1268 
1269 
1270  template <typename NumberType>
1271  Expression &
1272  Expression::operator*=(const NumberType &rhs)
1273  {
1274  *this = Expression(SymEngine::mul(get_RCP(), Expression(rhs).get_RCP()));
1275  return *this;
1276  }
1277 
1278 
1279  template <typename NumberType>
1280  Expression &
1281  Expression::operator/=(const NumberType &rhs)
1282  {
1283  *this = Expression(SymEngine::div(get_RCP(), Expression(rhs).get_RCP()));
1284  return *this;
1285  }
1286 
1287 
1288  template <typename ResultType>
1289  Expression::operator ResultType() const
1290  {
1291  return static_cast<ResultType>(get_expression());
1292  }
1293 
1294  } // namespace SD
1295 } // namespace Differentiation
1296 
1297 
1298 // Specializations from numbers.h
1299 // These are required in order to make the SD::Expression class a compatible
1300 // number for the Tensor and SymmetricTensor classes
1301 
1302 // Forward declarations:
1303 template <int rank_, int dim, typename Number>
1304 class Tensor;
1305 
1306 namespace internal
1307 {
1308  template <>
1309  struct NumberType<Differentiation::SD::Expression>
1310  {
1311  static const Differentiation::SD::Expression &
1312  value(const Differentiation::SD::Expression &t)
1313  {
1314  return t;
1315  }
1316 
1317  template <
1318  typename T,
1319  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1321  value(const T &t)
1322  {
1324  }
1325 
1326  template <
1327  typename T,
1328  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1330  value(T &&t)
1331  {
1333  }
1334 
1335  template <
1336  typename T,
1337  typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
1339  value(const std::complex<T> &t)
1340  {
1342  }
1343 
1344  template <typename T, int dim>
1346  value(const Tensor<0, dim, T> &t)
1347  {
1348  return Differentiation::SD::Expression(static_cast<T>(t));
1349  }
1350 
1351  template <typename T, int dim>
1353  value(const Tensor<0, dim, std::complex<T>> &t)
1354  {
1355  return Differentiation::SD::Expression(static_cast<std::complex<T>>(t));
1356  }
1357  };
1358 } // namespace internal
1359 
1360 
1361 namespace numbers
1362 {
1363  template <>
1364  inline bool
1366  {
1367  if (SymEngine::is_a_Number(value.get_value()))
1368  {
1369  const SymEngine::RCP<const SymEngine::Number> number_rcp =
1370  SymEngine::rcp_static_cast<const SymEngine::Number>(value.get_RCP());
1371  return number_rcp->is_zero();
1372  }
1373 
1374  return false;
1375  }
1376 
1377  template <>
1378  inline bool
1380  const Differentiation::SD::Expression &value_2)
1381  {
1382  return (value_1.get_value().__cmp__(value_2.get_value()) == 0);
1383  }
1384 
1385  template <>
1386  inline bool
1388  const Differentiation::SD::Expression &value_2)
1389  {
1390  return (value_1.get_value().__cmp__(value_2.get_value()) == -1);
1391  }
1392 } // namespace numbers
1393 
1394 
1395 # endif // DOXYGEN
1396 
1397 DEAL_II_NAMESPACE_CLOSE
1398 
1399 #endif // DEAL_II_WITH_SYMENGINE
1400 
1401 #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)
Expression operator>=(const Expression &lhs, const Expression &rhs)
Expression & parse(const std::string &expression)
bool value_is_zero(const Number &value)
Definition: numbers.h:876
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
bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:860
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:518
Expression operator<=(const Expression &lhs, const Expression &rhs)
#define Assert(cond, exc)
Definition: exceptions.h:1407
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
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
Expression operator*(Expression lhs, const Expression &rhs)
Expression operator!(const Expression &expression)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:884
Definition: mpi.h:90
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 &rhs)
Expression operator-(Expression lhs, const Expression &rhs)
Expression & operator+=(const Expression &rhs)
Expression operator!=(const Expression &lhs, const Expression &rhs)
Expression operator &(const Expression &lhs, const Expression &rhs)