Reference documentation for deal.II version GIT relicensing-767-g85e8698264 2024-05-26 14:00: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\}}\)
Loading...
Searching...
No Matches
symengine_number_types.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2019 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_differentiation_sd_symengine_number_types_h
16#define dealii_differentiation_sd_symengine_number_types_h
17
18#include <deal.II/base/config.h>
19
20#ifdef DEAL_II_WITH_SYMENGINE
21
22// Low level
23# include <symengine/basic.h>
24# include <symengine/dict.h>
25# include <symengine/symengine_exception.h>
26# include <symengine/symengine_rcp.h>
27
28// Number types
29# include <symengine/expression.h>
30# include <symengine/integer.h>
31# include <symengine/logic.h>
32# include <symengine/number.h>
33# include <symengine/rational.h>
34
35// Number operations
36# include <symengine/add.h>
37# include <symengine/functions.h>
38# include <symengine/mul.h>
39# include <symengine/pow.h>
40
41// Evaluation
42# include <symengine/eval.h>
43# include <symengine/eval_arb.h>
44# include <symengine/eval_double.h>
45
46// Differentiation
48# include <deal.II/base/numbers.h>
49
52
53# include <boost/serialization/split_member.hpp>
54
55# include <symengine/derivative.h>
56
57# include <algorithm>
58# include <memory>
59# include <sstream>
60# include <type_traits>
61# include <utility>
62# include <vector>
63
65
66namespace Differentiation
67{
68 namespace SD
69 {
81 std::string,
82 << "The string '" << arg1
83 << "' could not be parsed successfully. Are you sure that (1) it "
84 << "consists of legitimate operations and syntax, and (2) you've "
85 << "previously declared all symbolic variables that are present "
86 << "in the expression?");
87
177 {
178 public:
187 Expression();
188
196 explicit Expression(const bool value);
197
205 template <typename NumberType,
206 typename = std::enable_if_t<std::is_arithmetic_v<NumberType>>>
207 explicit Expression(const NumberType &value);
208
216 template <typename NumberType,
217 typename = std::enable_if_t<std::is_arithmetic_v<NumberType>>>
218 explicit Expression(const std::complex<NumberType> &value);
219
223 Expression(const SymEngine::integer_class &value);
224
231 template <typename NumberType,
232 typename = std::enable_if_t<std::is_integral_v<NumberType>>>
233 Expression(const NumberType &numerator, const NumberType &denominator);
234
238 Expression(const SymEngine::rational_class &value);
239
276 Expression(const Expression &condition,
277 const Expression &expression_if_true,
278 const Expression &expression_if_false);
279
301 Expression(const std::vector<std::pair<Expression, Expression>>
302 &condition_expression,
303 const Expression &expression_otherwise);
304
324 Expression(const std::vector<std::pair<Expression, Expression>>
325 &condition_expression);
326
327
334 Expression(const char *symbol);
335
348 Expression(const std::string &symb_expr,
349 const bool parse_as_expression = false);
350
357 Expression(const std::string &symbol_func,
358 const types::symbol_vector &arguments);
359
363 Expression(const Expression &rhs) = default;
364
372 explicit Expression(const SymEngine::Expression &rhs);
373
382 Expression(const SymEngine::RCP<const SymEngine::Basic> &rhs);
383
387 Expression(Expression &&rhs) = default;
388
397 Expression(SymEngine::RCP<const SymEngine::Basic> &&rhs);
398
402 virtual ~Expression() = default;
403
416 Expression &
417 parse(const std::string &expression);
418
425 std::ostream &
426 print(std::ostream &stream) const;
427
433 void
434 save(std::ostream &stream) const;
435
447 void
448 load(std::istream &stream);
449
458 template <class Archive>
459 void
460 save(Archive &archive, const unsigned int version) const;
461
474 template <class Archive>
475 void
476 load(Archive &archive, const unsigned int version);
477
478# ifdef DOXYGEN
493 template <class Archive>
494 void
495 serialize(Archive &archive, const unsigned int version);
496# else
497 // This macro defines the serialize() method that is compatible with
498 // the templated save() and load() method that have been implemented.
499 BOOST_SERIALIZATION_SPLIT_MEMBER()
500# endif
501
512 const SymEngine::Expression &
513 get_expression() const;
514
519 const SymEngine::Basic &
520 get_value() const;
521
526 const SymEngine::RCP<const SymEngine::Basic> &
527 get_RCP() const;
528
542 Expression &
543 operator=(const Expression &rhs);
544
551 Expression &
552 operator=(Expression &&rhs) noexcept;
553
560 Expression &
561 operator+=(const Expression &rhs);
562
569 Expression &
570 operator-=(const Expression &rhs);
571
578 Expression &
579 operator*=(const Expression &rhs);
580
587 Expression &
588 operator/=(const Expression &rhs);
589
603 template <typename NumberType>
604 Expression &
605 operator=(const NumberType &rhs);
606
616 operator-() const;
617
624 template <typename NumberType>
625 Expression &
626 operator+=(const NumberType &rhs);
627
634 template <typename NumberType>
635 Expression &
636 operator-=(const NumberType &rhs);
637
644 template <typename NumberType>
645 Expression &
646 operator*=(const NumberType &rhs);
647
654 template <typename NumberType>
655 Expression &
656 operator/=(const NumberType &rhs);
657
670 differentiate(const Expression &symbol) const;
671
678 const SymEngine::RCP<const SymEngine::Symbol> &symbol) const;
679
685 differentiate(const SymEngine::RCP<const SymEngine::Basic> &symbol) const;
686
707 substitute(const types::substitution_map &substitution_values) const;
708
726 substitute(const SymEngine::map_basic_basic &substitution_values) const;
727
738 substitute(const Expression &symbol, const Expression &value) const;
739
749 template <typename NumberType>
751 substitute(const Expression &symbol, const NumberType &value) const;
752
760 template <typename ReturnType>
761 ReturnType
763 const types::substitution_map &substitution_values) const;
764
776 template <typename ReturnType>
777 ReturnType
779 const SymEngine::map_basic_basic &substitution_values) const;
780
838 template <typename ResultType>
839 explicit operator ResultType() const;
840
845 explicit operator const SymEngine::Expression &() const;
846
851 operator const SymEngine::RCP<const SymEngine::Basic> &() const;
852
855 protected:
859 SymEngine::Expression &
861
862 private:
867 SymEngine::Expression expression;
868 };
869
880 template <>
881 struct is_symengine_number<Expression> : std::true_type
882 {};
883
884
890 template <>
891 struct is_sd_number<Expression> : std::true_type
892 {};
893
906 std::ostream &
907 operator<<(std::ostream &stream, const Expression &expression);
908
914 std::istream &
915 operator>>(std::istream &stream, Expression &expression);
916
930 operator==(const Expression &lhs, const Expression &rhs);
931
938 operator!=(const Expression &lhs, const Expression &rhs);
939
946 operator<(const Expression &lhs, const Expression &rhs);
947
954 operator>(const Expression &lhs, const Expression &rhs);
955
962 operator<=(const Expression &lhs, const Expression &rhs);
963
970 operator>=(const Expression &lhs, const Expression &rhs);
971
986 operator!(const Expression &expression);
987
995 operator&(const Expression &lhs, const Expression &rhs);
996
1004 operator|(const Expression &lhs, const Expression &rhs);
1005
1013 operator^(const Expression &lhs, const Expression &rhs);
1014
1023 operator&&(const Expression &lhs, const Expression &rhs);
1024
1033 operator||(const Expression &lhs, const Expression &rhs);
1034
1048 operator+(Expression lhs, const Expression &rhs);
1049
1056 operator-(Expression lhs, const Expression &rhs);
1057
1064 operator*(Expression lhs, const Expression &rhs);
1065
1072 operator/(Expression lhs, const Expression &rhs);
1073
1082 template <typename NumberType,
1083 typename = std::enable_if_t<
1084 std::is_constructible_v<Expression, NumberType>>>
1085 inline Expression
1086 operator+(const NumberType &lhs, const Expression &rhs)
1087 {
1088 return Expression(lhs) + rhs;
1089 }
1090
1099 template <typename NumberType,
1100 typename = std::enable_if_t<
1101 std::is_constructible_v<Expression, NumberType>>>
1102 inline Expression
1103 operator+(const Expression &lhs, const NumberType &rhs)
1104 {
1105 return lhs + Expression(rhs);
1106 }
1107
1116 template <typename NumberType,
1117 typename = std::enable_if_t<
1118 std::is_constructible_v<Expression, NumberType>>>
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 = std::enable_if_t<
1135 std::is_constructible_v<Expression, NumberType>>>
1136 inline Expression
1137 operator-(const Expression &lhs, const NumberType &rhs)
1138 {
1139 return lhs - Expression(rhs);
1140 }
1141
1150 template <typename NumberType,
1151 typename = std::enable_if_t<
1152 std::is_constructible_v<Expression, NumberType>>>
1153 inline Expression
1154 operator*(const NumberType &lhs, const Expression &rhs)
1155 {
1156 return Expression(lhs) * rhs;
1157 }
1158
1167 template <typename NumberType,
1168 typename = std::enable_if_t<
1169 std::is_constructible_v<Expression, NumberType>>>
1170 inline Expression
1171 operator*(const Expression &lhs, const NumberType &rhs)
1172 {
1173 return lhs * Expression(rhs);
1174 }
1175
1184 template <typename NumberType,
1185 typename = std::enable_if_t<
1186 std::is_constructible_v<Expression, NumberType>>>
1187 inline Expression
1188 operator/(const NumberType &lhs, const Expression &rhs)
1189 {
1190 return Expression(lhs) / rhs;
1191 }
1192
1201 template <typename NumberType,
1202 typename = std::enable_if_t<
1203 std::is_constructible_v<Expression, NumberType>>>
1204 inline Expression
1205 operator/(const Expression &lhs, const NumberType &rhs)
1206 {
1207 return lhs / Expression(rhs);
1208 }
1209
1212 } // namespace SD
1213} // namespace Differentiation
1214
1215
1216/* ---------------------- inline and template functions -------------------- */
1217
1218
1219# ifndef DOXYGEN
1220
1221
1222namespace Differentiation
1223{
1224 namespace SD
1225 {
1226 template <typename NumberType, typename>
1227 Expression::Expression(const NumberType &value)
1228 : expression(value)
1229 {}
1230
1231
1232 template <typename NumberType, typename>
1233 Expression::Expression(const std::complex<NumberType> &value)
1234 : expression(value)
1235 {}
1236
1237
1238 template <typename NumberType, typename>
1239 Expression::Expression(const NumberType &numerator,
1240 const NumberType &denominator)
1241 : expression(
1242 SymEngine::Rational::from_two_ints(*SymEngine::integer(numerator),
1243 *SymEngine::integer(denominator)))
1244 {}
1245
1246
1247 template <class Archive>
1248 void
1249 Expression::save(Archive &ar, const unsigned int /*version*/) const
1250 {
1251 std::stringstream sstream;
1252 sstream << *this;
1253 const std::string expr = sstream.str();
1254 ar &expr;
1255 }
1256
1257
1258 template <class Archive>
1259 void
1260 Expression::load(Archive &ar, const unsigned int /*version*/)
1261 {
1262 std::string expr;
1263 ar &expr;
1264 parse(expr);
1265 }
1266
1267
1268 template <typename NumberType>
1269 Expression
1270 Expression::substitute(const Expression &symbol,
1271 const NumberType &value) const
1272 {
1273 Assert(SymEngine::is_a<SymEngine::Symbol>(symbol.get_value()),
1274 ExcMessage(
1275 "Substitution with a number that does not represent a symbol."));
1276
1277 types::substitution_map sub_vals;
1278 sub_vals[symbol] = Expression(value);
1279 return substitute(sub_vals);
1280 }
1281
1282
1283 template <typename ReturnType>
1284 ReturnType
1285 Expression::substitute_and_evaluate(
1286 const types::substitution_map &substitution_values) const
1287 {
1288 return static_cast<ReturnType>(substitute(substitution_values));
1289 }
1290
1291
1292 template <typename ReturnType>
1293 ReturnType
1294 Expression::substitute_and_evaluate(
1295 const SymEngine::map_basic_basic &substitution_values) const
1296 {
1297 return static_cast<ReturnType>(substitute(substitution_values));
1298 }
1299
1300
1301 template <typename NumberType>
1302 Expression &
1303 Expression::operator=(const NumberType &rhs)
1304 {
1305 *this = Expression(rhs);
1306 return *this;
1307 }
1308
1309
1310 template <typename NumberType>
1311 Expression &
1312 Expression::operator+=(const NumberType &rhs)
1313 {
1314 *this = Expression(SymEngine::add(get_RCP(), Expression(rhs).get_RCP()));
1315 return *this;
1316 }
1317
1318
1319 template <typename NumberType>
1320 Expression &
1321 Expression::operator-=(const NumberType &rhs)
1322 {
1323 *this = Expression(SymEngine::sub(get_RCP(), Expression(rhs).get_RCP()));
1324 return *this;
1325 }
1326
1327
1328 template <typename NumberType>
1329 Expression &
1330 Expression::operator*=(const NumberType &rhs)
1331 {
1332 *this = Expression(SymEngine::mul(get_RCP(), Expression(rhs).get_RCP()));
1333 return *this;
1334 }
1335
1336
1337 template <typename NumberType>
1338 Expression &
1339 Expression::operator/=(const NumberType &rhs)
1340 {
1341 *this = Expression(SymEngine::div(get_RCP(), Expression(rhs).get_RCP()));
1342 return *this;
1343 }
1344
1345
1346 template <typename ResultType>
1347 Expression::operator ResultType() const
1348 {
1349 return static_cast<ResultType>(get_expression());
1350 }
1351
1352 } // namespace SD
1353} // namespace Differentiation
1354
1355
1356// Specializations from numbers.h
1357// These are required in order to make the SD::Expression class a compatible
1358// number for the Tensor and SymmetricTensor classes
1359
1360// Forward declarations:
1361template <int rank_, int dim, typename Number>
1362class Tensor;
1363
1364namespace internal
1365{
1366 template <>
1367 struct NumberType<Differentiation::SD::Expression>
1368 {
1369 static const Differentiation::SD::Expression &
1371 {
1372 return t;
1373 }
1374
1375 template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
1377 value(const T &t)
1378 {
1380 }
1381
1382 template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
1384 value(T &&t)
1385 {
1387 }
1388
1389 template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
1391 value(const std::complex<T> &t)
1392 {
1394 }
1395
1396 template <typename T, int dim>
1398 value(const Tensor<0, dim, T> &t)
1399 {
1400 return Differentiation::SD::Expression(static_cast<T>(t));
1401 }
1402
1403 template <typename T, int dim>
1405 value(const Tensor<0, dim, std::complex<T>> &t)
1406 {
1407 return Differentiation::SD::Expression(static_cast<std::complex<T>>(t));
1408 }
1409 };
1410} // namespace internal
1411
1412
1413namespace numbers
1414{
1415 template <>
1416 inline bool
1418 {
1419 if (SymEngine::is_a_Number(value.get_value()))
1420 {
1421 const SymEngine::RCP<const SymEngine::Number> number_rcp =
1422 SymEngine::rcp_static_cast<const SymEngine::Number>(value.get_RCP());
1423 return number_rcp->is_zero();
1424 }
1425
1426 return false;
1427 }
1428
1429 template <>
1430 inline bool
1432 const Differentiation::SD::Expression &value_2)
1433 {
1434 return (value_1.get_value().__cmp__(value_2.get_value()) == 0);
1435 }
1436
1437 template <>
1438 inline bool
1440 const Differentiation::SD::Expression &value_2)
1441 {
1442 return (value_1.get_value().__cmp__(value_2.get_value()) == -1);
1443 }
1444} // namespace numbers
1445
1446
1447# endif // DOXYGEN
1448
1450
1451#endif // DEAL_II_WITH_SYMENGINE
1452
1453#endif // dealii_differentiation_sd_symengine_number_types_h
Expression & operator/=(const Expression &rhs)
Expression & parse(const std::string &expression)
Expression & operator-=(const NumberType &rhs)
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 & operator/=(const NumberType &rhs)
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)
std::ostream & print(std::ostream &stream) const
Expression substitute(const types::substitution_map &substitution_values) const
Expression & operator=(const NumberType &rhs)
void save(Archive &archive, const unsigned int version) const
Expression(const Expression &rhs)=default
void save(std::ostream &stream) const
Expression & operator*=(const NumberType &rhs)
ReturnType substitute_and_evaluate(const SymEngine::map_basic_basic &substitution_values) const
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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcSymEngineParserError(std::string arg1)
#define Assert(cond, exc)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
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)
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)
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:936
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:920
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:944