Reference documentation for deal.II version GIT 00e6fe71c9 2023-06-04 19:35:01+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 - 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
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 <
207  typename NumberType,
208  typename = std::enable_if_t<std::is_arithmetic<NumberType>::value>>
209  explicit Expression(const NumberType &value);
210 
218  template <
219  typename NumberType,
220  typename = std::enable_if_t<std::is_arithmetic<NumberType>::value>>
221  explicit Expression(const std::complex<NumberType> &value);
222 
226  Expression(const SymEngine::integer_class &value);
227 
234  template <
235  typename NumberType,
236  typename = std::enable_if_t<std::is_integral<NumberType>::value>>
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 
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 
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 
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 
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 
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 
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 
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 
859  protected:
863  SymEngine::Expression &
864  get_expression();
865 
866  private:
871  SymEngine::Expression expression;
872  };
873 
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 
910  std::ostream &
911  operator<<(std::ostream &stream, const Expression &expression);
912 
918  std::istream &
919  operator>>(std::istream &stream, Expression &expression);
920 
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 
989  Expression
990  operator!(const Expression &expression);
991 
998  Expression
999  operator&(const Expression &lhs, const Expression &rhs);
1000 
1007  Expression
1008  operator|(const Expression &lhs, const Expression &rhs);
1009 
1016  Expression
1017  operator^(const Expression &lhs, const Expression &rhs);
1018 
1026  Expression
1027  operator&&(const Expression &lhs, const Expression &rhs);
1028 
1036  Expression
1037  operator||(const Expression &lhs, const Expression &rhs);
1038 
1051  Expression
1052  operator+(Expression lhs, const Expression &rhs);
1053 
1059  Expression
1060  operator-(Expression lhs, const Expression &rhs);
1061 
1067  Expression
1068  operator*(Expression lhs, const Expression &rhs);
1069 
1075  Expression
1076  operator/(Expression lhs, const Expression &rhs);
1077 
1086  template <typename NumberType,
1087  typename = std::enable_if_t<
1088  std::is_constructible<Expression, NumberType>::value>>
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 = std::enable_if_t<
1105  std::is_constructible<Expression, NumberType>::value>>
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 = std::enable_if_t<
1122  std::is_constructible<Expression, NumberType>::value>>
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 = std::enable_if_t<
1139  std::is_constructible<Expression, NumberType>::value>>
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 = std::enable_if_t<
1156  std::is_constructible<Expression, NumberType>::value>>
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 = std::enable_if_t<
1173  std::is_constructible<Expression, NumberType>::value>>
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 = std::enable_if_t<
1190  std::is_constructible<Expression, NumberType>::value>>
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 = std::enable_if_t<
1207  std::is_constructible<Expression, NumberType>::value>>
1208  inline Expression
1209  operator/(const Expression &lhs, const NumberType &rhs)
1210  {
1211  return lhs / Expression(rhs);
1212  }
1213 
1216  } // namespace SD
1217 } // namespace Differentiation
1218 
1219 
1220 /* ---------------------- inline and template functions -------------------- */
1221 
1222 
1223 # ifndef DOXYGEN
1224 
1225 
1226 namespace 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
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
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:
1365 template <int rank_, int dim, typename Number>
1366 class Tensor;
1367 
1368 namespace 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 <typename T,
1380  typename = std::enable_if_t<std::is_arithmetic<T>::value>>
1382  value(const T &t)
1383  {
1385  }
1386 
1387  template <typename T,
1388  typename = std::enable_if_t<std::is_arithmetic<T>::value>>
1390  value(T &&t)
1391  {
1393  }
1394 
1395  template <typename T,
1396  typename = std::enable_if_t<std::is_arithmetic<T>::value>>
1398  value(const std::complex<T> &t)
1399  {
1401  }
1402 
1403  template <typename T, int dim>
1405  value(const Tensor<0, dim, T> &t)
1406  {
1407  return Differentiation::SD::Expression(static_cast<T>(t));
1408  }
1409 
1410  template <typename T, int dim>
1412  value(const Tensor<0, dim, std::complex<T>> &t)
1413  {
1414  return Differentiation::SD::Expression(static_cast<std::complex<T>>(t));
1415  }
1416  };
1417 } // namespace internal
1418 
1419 
1420 namespace numbers
1421 {
1422  template <>
1423  inline bool
1425  {
1426  if (SymEngine::is_a_Number(value.get_value()))
1427  {
1428  const SymEngine::RCP<const SymEngine::Number> number_rcp =
1429  SymEngine::rcp_static_cast<const SymEngine::Number>(value.get_RCP());
1430  return number_rcp->is_zero();
1431  }
1432 
1433  return false;
1434  }
1435 
1436  template <>
1437  inline bool
1439  const Differentiation::SD::Expression &value_2)
1440  {
1441  return (value_1.get_value().__cmp__(value_2.get_value()) == 0);
1442  }
1443 
1444  template <>
1445  inline bool
1447  const Differentiation::SD::Expression &value_2)
1448  {
1449  return (value_1.get_value().__cmp__(value_2.get_value()) == -1);
1450  }
1451 } // namespace numbers
1452 
1453 
1454 # endif // DOXYGEN
1455 
1457 
1458 #endif // DEAL_II_WITH_SYMENGINE
1459 
1460 #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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcSymEngineParserError(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1614
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
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:939
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:923
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:947
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:702