Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50: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\}}\)
symengine_number_types.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 - 2022 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 
178  {
179  public:
188  Expression();
189 
197  explicit Expression(const bool value);
198 
206  template <typename NumberType,
207  typename = std::enable_if_t<std::is_arithmetic_v<NumberType>>>
208  explicit Expression(const NumberType &value);
209 
217  template <typename NumberType,
218  typename = std::enable_if_t<std::is_arithmetic_v<NumberType>>>
219  explicit Expression(const std::complex<NumberType> &value);
220 
224  Expression(const SymEngine::integer_class &value);
225 
232  template <typename NumberType,
233  typename = std::enable_if_t<std::is_integral_v<NumberType>>>
234  Expression(const NumberType &numerator, const NumberType &denominator);
235 
239  Expression(const SymEngine::rational_class &value);
240 
277  Expression(const Expression &condition,
278  const Expression &expression_if_true,
279  const Expression &expression_if_false);
280 
302  Expression(const std::vector<std::pair<Expression, Expression>>
303  &condition_expression,
304  const Expression &expression_otherwise);
305 
325  Expression(const std::vector<std::pair<Expression, Expression>>
326  &condition_expression);
327 
328 
335  Expression(const char *symbol);
336 
349  Expression(const std::string &symb_expr,
350  const bool parse_as_expression = false);
351 
358  Expression(const std::string &symbol_func,
359  const types::symbol_vector &arguments);
360 
364  Expression(const Expression &rhs) = default;
365 
373  explicit Expression(const SymEngine::Expression &rhs);
374 
383  Expression(const SymEngine::RCP<const SymEngine::Basic> &rhs);
384 
388  Expression(Expression &&rhs) = default;
389 
398  Expression(SymEngine::RCP<const SymEngine::Basic> &&rhs);
399 
403  virtual ~Expression() = default;
404 
417  Expression &
418  parse(const std::string &expression);
419 
426  std::ostream &
427  print(std::ostream &stream) const;
428 
434  void
435  save(std::ostream &stream) const;
436 
448  void
449  load(std::istream &stream);
450 
459  template <class Archive>
460  void
461  save(Archive &archive, const unsigned int version) const;
462 
475  template <class Archive>
476  void
477  load(Archive &archive, const unsigned int version);
478 
479 # ifdef DOXYGEN
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 
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 
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 
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 
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 
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 
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 
856  protected:
860  SymEngine::Expression &
861  get_expression();
862 
863  private:
868  SymEngine::Expression expression;
869  };
870 
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 
907  std::ostream &
908  operator<<(std::ostream &stream, const Expression &expression);
909 
915  std::istream &
916  operator>>(std::istream &stream, Expression &expression);
917 
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 
986  Expression
987  operator!(const Expression &expression);
988 
995  Expression
996  operator&(const Expression &lhs, const Expression &rhs);
997 
1004  Expression
1005  operator|(const Expression &lhs, const Expression &rhs);
1006 
1013  Expression
1014  operator^(const Expression &lhs, const Expression &rhs);
1015 
1023  Expression
1024  operator&&(const Expression &lhs, const Expression &rhs);
1025 
1033  Expression
1034  operator||(const Expression &lhs, const Expression &rhs);
1035 
1048  Expression
1049  operator+(Expression lhs, const Expression &rhs);
1050 
1056  Expression
1057  operator-(Expression lhs, const Expression &rhs);
1058 
1064  Expression
1065  operator*(Expression lhs, const Expression &rhs);
1066 
1072  Expression
1073  operator/(Expression lhs, const Expression &rhs);
1074 
1083  template <typename NumberType,
1084  typename = std::enable_if_t<
1085  std::is_constructible_v<Expression, NumberType>>>
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 = std::enable_if_t<
1102  std::is_constructible_v<Expression, NumberType>>>
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 = std::enable_if_t<
1119  std::is_constructible_v<Expression, NumberType>>>
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 = std::enable_if_t<
1136  std::is_constructible_v<Expression, NumberType>>>
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 = std::enable_if_t<
1153  std::is_constructible_v<Expression, NumberType>>>
1154  inline Expression
1155  operator*(const NumberType &lhs, const Expression &rhs)
1156  {
1157  return Expression(lhs) * rhs;
1158  }
1159 
1168  template <typename NumberType,
1169  typename = std::enable_if_t<
1170  std::is_constructible_v<Expression, NumberType>>>
1171  inline Expression
1172  operator*(const Expression &lhs, const NumberType &rhs)
1173  {
1174  return lhs * Expression(rhs);
1175  }
1176 
1185  template <typename NumberType,
1186  typename = std::enable_if_t<
1187  std::is_constructible_v<Expression, NumberType>>>
1188  inline Expression
1189  operator/(const NumberType &lhs, const Expression &rhs)
1190  {
1191  return Expression(lhs) / rhs;
1192  }
1193 
1202  template <typename NumberType,
1203  typename = std::enable_if_t<
1204  std::is_constructible_v<Expression, NumberType>>>
1205  inline Expression
1206  operator/(const Expression &lhs, const NumberType &rhs)
1207  {
1208  return lhs / Expression(rhs);
1209  }
1210 
1213  } // namespace SD
1214 } // namespace Differentiation
1215 
1216 
1217 /* ---------------------- inline and template functions -------------------- */
1218 
1219 
1220 # ifndef DOXYGEN
1221 
1222 
1223 namespace Differentiation
1224 {
1225  namespace SD
1226  {
1227  template <typename NumberType, typename>
1228  Expression::Expression(const NumberType &value)
1229  : expression(value)
1230  {}
1231 
1232 
1233  template <typename NumberType, typename>
1234  Expression::Expression(const std::complex<NumberType> &value)
1235  : expression(value)
1236  {}
1237 
1238 
1239  template <typename NumberType, typename>
1240  Expression::Expression(const NumberType &numerator,
1241  const NumberType &denominator)
1242  : expression(
1243  SymEngine::Rational::from_two_ints(*SymEngine::integer(numerator),
1244  *SymEngine::integer(denominator)))
1245  {}
1246 
1247 
1248  template <class Archive>
1249  void
1250  Expression::save(Archive &ar, const unsigned int /*version*/) const
1251  {
1252  std::stringstream sstream;
1253  sstream << *this;
1254  const std::string expr = sstream.str();
1255  ar &expr;
1256  }
1257 
1258 
1259  template <class Archive>
1260  void
1261  Expression::load(Archive &ar, const unsigned int /*version*/)
1262  {
1263  std::string expr;
1264  ar &expr;
1265  parse(expr);
1266  }
1267 
1268 
1269  template <typename NumberType>
1270  Expression
1271  Expression::substitute(const Expression &symbol,
1272  const NumberType &value) const
1273  {
1274  Assert(SymEngine::is_a<SymEngine::Symbol>(symbol.get_value()),
1275  ExcMessage(
1276  "Substitution with a number that does not represent a symbol."));
1277 
1278  types::substitution_map sub_vals;
1279  sub_vals[symbol] = Expression(value);
1280  return substitute(sub_vals);
1281  }
1282 
1283 
1284  template <typename ReturnType>
1285  ReturnType
1287  const types::substitution_map &substitution_values) const
1288  {
1289  return static_cast<ReturnType>(substitute(substitution_values));
1290  }
1291 
1292 
1293  template <typename ReturnType>
1294  ReturnType
1296  const SymEngine::map_basic_basic &substitution_values) const
1297  {
1298  return static_cast<ReturnType>(substitute(substitution_values));
1299  }
1300 
1301 
1302  template <typename NumberType>
1303  Expression &
1304  Expression::operator=(const NumberType &rhs)
1305  {
1306  *this = Expression(rhs);
1307  return *this;
1308  }
1309 
1310 
1311  template <typename NumberType>
1312  Expression &
1313  Expression::operator+=(const NumberType &rhs)
1314  {
1315  *this = Expression(SymEngine::add(get_RCP(), Expression(rhs).get_RCP()));
1316  return *this;
1317  }
1318 
1319 
1320  template <typename NumberType>
1321  Expression &
1322  Expression::operator-=(const NumberType &rhs)
1323  {
1324  *this = Expression(SymEngine::sub(get_RCP(), Expression(rhs).get_RCP()));
1325  return *this;
1326  }
1327 
1328 
1329  template <typename NumberType>
1330  Expression &
1331  Expression::operator*=(const NumberType &rhs)
1332  {
1333  *this = Expression(SymEngine::mul(get_RCP(), Expression(rhs).get_RCP()));
1334  return *this;
1335  }
1336 
1337 
1338  template <typename NumberType>
1339  Expression &
1340  Expression::operator/=(const NumberType &rhs)
1341  {
1342  *this = Expression(SymEngine::div(get_RCP(), Expression(rhs).get_RCP()));
1343  return *this;
1344  }
1345 
1346 
1347  template <typename ResultType>
1348  Expression::operator ResultType() const
1349  {
1350  return static_cast<ResultType>(get_expression());
1351  }
1352 
1353  } // namespace SD
1354 } // namespace Differentiation
1355 
1356 
1357 // Specializations from numbers.h
1358 // These are required in order to make the SD::Expression class a compatible
1359 // number for the Tensor and SymmetricTensor classes
1360 
1361 // Forward declarations:
1362 template <int rank_, int dim, typename Number>
1363 class Tensor;
1364 
1365 namespace internal
1366 {
1367  template <>
1368  struct NumberType<Differentiation::SD::Expression>
1369  {
1370  static const Differentiation::SD::Expression &
1372  {
1373  return t;
1374  }
1375 
1376  template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
1378  value(const T &t)
1379  {
1381  }
1382 
1383  template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
1385  value(T &&t)
1386  {
1388  }
1389 
1390  template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
1392  value(const std::complex<T> &t)
1393  {
1395  }
1396 
1397  template <typename T, int dim>
1399  value(const Tensor<0, dim, T> &t)
1400  {
1401  return Differentiation::SD::Expression(static_cast<T>(t));
1402  }
1403 
1404  template <typename T, int dim>
1406  value(const Tensor<0, dim, std::complex<T>> &t)
1407  {
1408  return Differentiation::SD::Expression(static_cast<std::complex<T>>(t));
1409  }
1410  };
1411 } // namespace internal
1412 
1413 
1414 namespace numbers
1415 {
1416  template <>
1417  inline bool
1419  {
1420  if (SymEngine::is_a_Number(value.get_value()))
1421  {
1422  const SymEngine::RCP<const SymEngine::Number> number_rcp =
1423  SymEngine::rcp_static_cast<const SymEngine::Number>(value.get_RCP());
1424  return number_rcp->is_zero();
1425  }
1426 
1427  return false;
1428  }
1429 
1430  template <>
1431  inline bool
1433  const Differentiation::SD::Expression &value_2)
1434  {
1435  return (value_1.get_value().__cmp__(value_2.get_value()) == 0);
1436  }
1437 
1438  template <>
1439  inline bool
1441  const Differentiation::SD::Expression &value_2)
1442  {
1443  return (value_1.get_value().__cmp__(value_2.get_value()) == -1);
1444  }
1445 } // namespace numbers
1446 
1447 
1448 # endif // DOXYGEN
1449 
1451 
1452 #endif // DEAL_II_WITH_SYMENGINE
1453 
1454 #endif // dealii_differentiation_sd_symengine_number_types_h
Expression & operator/=(const Expression &rhs)
Expression & parse(const std::string &expression)
Expression substitute(const Expression &symbol, const NumberType &value) const
Expression & operator=(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(const NumberType &numerator, const NumberType &denominator)
Expression(Expression &&rhs)=default
Expression & operator*=(const Expression &rhs)
void serialize(Archive &archive, const unsigned int version)
Expression & operator-=(const NumberType &rhs)
Expression & operator+=(const NumberType &rhs)
std::ostream & print(std::ostream &stream) const
Expression substitute(const types::substitution_map &substitution_values) const
void save(Archive &archive, const unsigned int version) const
Expression & operator*=(const NumberType &rhs)
Expression(const Expression &rhs)=default
void save(std::ostream &stream) const
ReturnType substitute_and_evaluate(const SymEngine::map_basic_basic &substitution_values) const
Expression & operator/=(const NumberType &rhs)
Expression & operator=(const NumberType &rhs)
Expression & operator-=(const Expression &rhs)
const SymEngine::Basic & get_value() const
const SymEngine::Expression & get_expression() const
Expression & operator+=(const Expression &rhs)
void load(Archive &archive, const unsigned int version)
Expression(const std::complex< NumberType > &value)
Expression differentiate(const Expression &symbol) const
Definition: tensor.h:516
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcSymEngineParserError(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:517
static ::ExceptionBase & ExcMessage(std::string arg1)
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)
ValueType substitute_and_evaluate(const Expression &expression, const types::substitution_map &substitution_map)
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)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
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)
static const char T
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:937
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:921
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:945
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:703