--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_differentiation_sd_symengine_number_types_h
+#define dealii_differentiation_sd_symengine_number_types_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SYMENGINE
+
+// Low level
+# include <symengine/basic.h>
+# include <symengine/dict.h>
+# include <symengine/symengine_exception.h>
+# include <symengine/symengine_rcp.h>
+
+// Number types
+# include <symengine/expression.h>
+# include <symengine/integer.h>
+# include <symengine/logic.h>
+# include <symengine/rational.h>
+
+// Number operations
+# include <symengine/add.h>
+# include <symengine/functions.h>
+# include <symengine/mul.h>
+# include <symengine/pow.h>
+
+// Evaluation
+# include <symengine/eval.h>
+# include <symengine/eval_arb.h>
+# include <symengine/eval_double.h>
+
+// Differentiation
+# include <deal.II/base/exceptions.h>
+
+# include <deal.II/differentiation/sd/symengine_number_traits.h>
+# include <deal.II/differentiation/sd/symengine_types.h>
+
+# include <symengine/derivative.h>
+
+# include <algorithm>
+# include <memory>
+# include <sstream>
+# include <type_traits>
+# include <utility>
+# include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Differentiation
+{
+ namespace SD
+ {
+ namespace SE = ::SymEngine;
+
+ /**
+ * @addtogroup Exceptions
+ * @{
+ */
+
+ /**
+ * An exception to indicate that the string sent to the SymEngine
+ * parser is not valid.
+ */
+ DeclException1(
+ ExcSymEngineParserError,
+ std::string,
+ << "The string '" << arg1
+ << "' could not be parsed successfully. Are you sure that (1) it "
+ << "consists of legitimate operations and syntax, and (2) you've "
+ << "previously declared all symbolic variables that are present "
+ << "in the expression?");
+
+ //@}
+
+ /**
+ * A class to wrap SymEngine expressions.
+ *
+ * With this number type, SymEngine numbers can be used to perform scalar
+ * and tensor mathematics in deal.II. It (or the SymEngine::Expression
+ * class, of which it stores an instance) therefore forms the basis of
+ * symbolic computation via SymEngine in deal.II. With it one can
+ * perform symbolic differentiation and subsequent substitution with both
+ * scalars and deal.II's native Tensor and SymmetricTensor types.
+ *
+ * The symbolic features that this class supports includes:
+ * - expression parsing,
+ * - comparison operations,
+ * - logical operations,
+ * - math operations,
+ * - conditional expression construction,
+ * - differentiation,
+ * - substitution (partial and complete), and
+ * - serialization.
+ *
+ * A simple example of how this class may be used is as follows:
+ * @code
+ * // Constructing a symbolic expression:
+ * // This is a symbol, which we will treat as an argument to a symbolic
+ * // function.
+ * const Expression x("x");
+ * const Expression y("y");
+ * // This is a symbolic expression, which is an expression constructed
+ * // from individual symbols.
+ * const Expression f = (x + y)*(x + y);
+ *
+ * // Value substitution
+ * types::substitution_map substitution_map;
+ * substitution_map[x] = Expression(1);
+ * substitution_map[y] = Expression(2.5);
+ * const double evaluated_f =
+ * f.substitute_and_evaluate<double>(substitution_map);
+ * // We could also have performed substitution of each individual
+ * // argument, if we wanted to. This means that one can partially
+ * // substitute an expression at any time.
+ * @endcode
+ *
+ * A more intricate example of conditional evaluation is as follows:
+ * @code
+ * // Construct symbolic expressions
+ * const Expression x("x");
+ * const Expression y("y");
+ * const Expression f_plus = (x + y)*(x + y);
+ * // Parsing expressions from a string is also possible. Its arguments
+ * // must have been previously declared through (and in scope).
+ * const Expression f_minus ("(x-y)*(x-y)", true);
+ *
+ * // Constructing a conditional expression
+ * const SD_number_t f((x > Expression(0.0)), f_plus, f_minus);
+ *
+ * // Value substitution
+ * types::substitution_map substitution_map;
+ * substitution_map[x] = Expression(1);
+ * substitution_map[y] = Expression(2.5);
+ * const double evaluated_f =
+ * f.substitute_and_evaluate<double>(substitution_map);
+ * // Since the substituted value for x was greater than zero, we expect
+ * // that the returned result now in evaluated_f was evaluated from
+ * // the function f_plus.
+ * @endcode
+ *
+ * Lastly, here is an example using symbolic differentiation:
+ * @code
+ * // Construct symbolic expressions
+ * const Expression x("x");
+ * const Expression f("x**2", true);
+ *
+ * // Now perform differentiation. Specifically, we differentiate the
+ * // function "f" with respect to the symbolic variable "x".
+ * // The result should be the expression "2*x".
+ * const Expression df_dx = f.differentiate(x);
+ *
+ * // Value substitution
+ * types::substitution_map substitution_map;
+ * substitution_map[x] = Expression(10.0);
+ * const double evaluated_df_dx =
+ * evaluated_df_dx.substitute_and_evaluate<double>(substitution_map);
+ * // We can expect the above to evaluate to "2*10" which is,
+ * // of course, the numeric value 20.
+ * @endcode
+ *
+ * @author Jean-Paul Pelteret, 2019
+ */
+ class Expression
+ {
+ public:
+ /**
+ * @name Constructors
+ */
+ //@{
+
+ /**
+ * Default constructor.
+ */
+ Expression();
+
+ /**
+ * Constructor for boolean types.
+ *
+ * @note This constructor is marked as explicit so that there are no
+ * potential ambiguities related to implicit conversions in either user
+ * code or math functions that are loaded into the standard namespace.
+ */
+ explicit Expression(const bool &value);
+
+ /**
+ * Constructor for arithmetic number types.
+ *
+ * @note This constructor is marked as explicit so that there are no
+ * potential ambiguities related to implicit conversions in either user
+ * code or math functions that are loaded into the standard namespace.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_arithmetic<NumberType>::value>::type>
+ explicit Expression(const NumberType &value);
+
+ /**
+ * Constructor for complex numbers templated on arithmetic number types.
+ *
+ * @note This constructor is marked as explicit so that there are no
+ * potential ambiguities related to implicit conversions in either user
+ * code or math functions that are loaded into the standard namespace.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_arithmetic<NumberType>::value>::type>
+ explicit Expression(const std::complex<NumberType> &value);
+
+ /**
+ * Constructor for integer types.
+ */
+ Expression(const SE::integer_class &value);
+
+ /**
+ * Constructor for rational types.
+ *
+ * It is expected that both the @p numerator and @p denominator
+ * be integral types.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_integral<NumberType>::value>::type>
+ Expression(const NumberType &numerator, const NumberType &denominator);
+
+ /**
+ * Constructor for rational types.
+ */
+ Expression(const SE::rational_class &value);
+
+ /**
+ * Constructor for a piecewise defined function.
+ *
+ * The generated expression may be interpreted as the result of the
+ * teniary operator, i.e.
+ * <code>(condition ? expression_if_true : expression_if_false)</code>.
+ *
+ * The @p condition can be any expression that renders an expression that
+ * is convertible to a SymEngine::Boolean operator. This includes:
+ * - the logical operators of the deal.II Expression class
+ * - SymEngine::boolean()
+ * - SymEngine::contains()
+ * - SymEngine::Eq()
+ * - SymEngine::Ne()
+ * - SymEngine::Ge()
+ * - SymEngine::Gt()
+ * - SymEngine::Le()
+ * - SymEngine::Lt()
+ * - SymEngine::logical_and()
+ * - SymEngine::logical_nand()
+ * - SymEngine::logical_or()
+ * - SymEngine::logical_not()
+ * - SymEngine::logical_nor()
+ * - SymEngine::logical_xor()
+ * - SymEngine::logical_xnor()
+ * - ...
+ *
+ * An example of this constructor's use is as follows:
+ * @code
+ * const Expression x("x");
+ * const Expression y("y");
+ *
+ * // Construct a conditional expression using the symbolic variables.
+ * const Expression f ((x < Expression(0.0)), x+y, x-y);
+ * @endcode
+ */
+ Expression(const Expression &condition,
+ const Expression &expression_if_true,
+ const Expression &expression_if_false);
+
+ /**
+ * Constructor for a piecewise defined function.
+ *
+ * The generated expression may be interpreted as the result of the
+ * set of nested if-elseif-else statements, i.e. (in pseudo-code)
+ * @code
+ * if (condition_expression[0].first == true)
+ * return condition_expression[0].second;
+ * else if (condition_expression[1].first == true)
+ * return condition_expression[1].second;
+ * else if (...)
+ * return ...;
+ * else
+ * return expression_otherwise;
+ * @endcode
+ * if the input vector has more than 2 elements.
+ *
+ * This variant takes the piecewise evaluated conditions and its results
+ * as the first argument, and the default return value as the second
+ * argument.
+ */
+ Expression(const std::vector<std::pair<Expression, Expression>>
+ & condition_expression,
+ const Expression &expression_otherwise);
+
+ /**
+ * Constructor for a piecewise defined function.
+ *
+ * The generated expression may be interpreted as the result of the
+ * set of nested if-elseif statements, i.e. (in pseudo-code)
+ * @code
+ * if (condition_expression[0].first == true)
+ * return condition_expression[0].second;
+ * else if (condition_expression[1].first == true)
+ * return condition_expression[1].second;
+ * else if (...)
+ * return ...;
+ * @endcode
+ * if the input vector has more than 2 elements.
+ *
+ * This variant takes only the piecewise evaluated conditions and its
+ * results. If none of the conditions are met upon evaluation then the
+ * returned result will be NaN.
+ */
+ Expression(const std::vector<std::pair<Expression, Expression>>
+ &condition_expression);
+
+
+ /**
+ * Constructor for symbolic types.
+ *
+ * This constructor initializes a symbolic type with a character array
+ * representing its symbolic value.
+ */
+ Expression(const char *symbol);
+
+ /**
+ * Constructor for symbolic types.
+ *
+ * This constructor initializes a symbolic type with a string
+ * representing its symbolic value.
+ * If the @p parse_as_expression flag is <tt>false</tt>, then the
+ * @p symb_expr (potentially composed of mulitple characters) will be
+ * interpreted as a single symbol.
+ * If the @p parse_as_expression flag is <tt>true</tt>, then the
+ * @p symb_expr will be parsed as a symbolic expression (potentially
+ * composed of multiple symbols, constants, etc.).
+ */
+ Expression(const std::string &symb_expr,
+ const bool parse_as_expression = false);
+
+ /**
+ * Constructor for function symbol types.
+ *
+ * This constructor initializes a function symbol with a string
+ * representing its symbolic name.
+ */
+ Expression(const std::string & symbol_func,
+ const types::symbol_vector &arguments);
+
+ /**
+ * Copy constructor.
+ */
+ Expression(const Expression &rhs) = default;
+
+ /**
+ * Copy constructor.
+ *
+ * @note This constructor is marked as explicit to prevent any ambiguities
+ * from implicit conversion when both the deal.II and SymEngine namespaces
+ * are imported.
+ */
+ explicit Expression(const SE::Expression &rhs);
+
+ /**
+ * Copy constructor.
+ *
+ * This allows us to create our class straight from the results of
+ * SymEngine operations. This is especially important for operations like
+ * "diff", because the returned result is not primitive, but rather a set
+ * of compound operations.
+ */
+ Expression(const SE::RCP<const SE::Basic> &rhs);
+
+ /**
+ * Move constructor.
+ */
+ Expression(Expression &&rhs) = default;
+
+ /**
+ * Move constructor.
+ *
+ * This allows us to create our class straight from the results of
+ * SymEngine operations. This is especially important for operations like
+ * "diff", because the returned result is not primitive, but rather a set
+ * of compound operations.
+ */
+ Expression(SE::RCP<const SE::Basic> &&rhs);
+
+ /**
+ * Destructor.
+ */
+ virtual ~Expression() = default;
+
+ //@}
+
+ /**
+ * Utilities
+ */
+ //@{
+
+ /**
+ * Parse an expression from a string representing a symbolic @p expression.
+ * This overwrites any existing value or expression that this object
+ * represents.
+ */
+ Expression &
+ parse(const std::string &expression);
+
+ /**
+ * Print the value stored by this object.
+ *
+ * Since the stored value could be one of a number of types, we leave
+ * SymEngine to cast and output the correct representation of the data.
+ */
+ std::ostream &
+ print(std::ostream &stream) const;
+
+ /**
+ * Save the value stored by this object to the @p stream.
+ *
+ * Each expression will be saved on a new line of the @p stream.
+ */
+ void
+ save(std::ostream &stream) const;
+
+ /**
+ * Load the value stored in the @p stream into this object.
+ *
+ * It is expected that each expression appears on its own on single
+ * line of @p stream.
+ *
+ * @note When loading a symbolic expression, it is imperative that
+ * you first create or load all of the symbolic variables used in
+ * the saved expression.
+ */
+ void
+ load(std::istream &stream);
+
+ /**
+ * Write and read the data of this object from a stream for the purpose
+ * of serialization.
+ *
+ * This effecively saves or loads the value stored into/out of the
+ * @p archive with the given @p version number into this object.
+ * If deserializing data, then the previous contents of this object
+ * are thrown away.
+ *
+ * @note When deserializing a symbolic expression, it is imperative that
+ * you first create or deserialize all of the symbolic variables used in
+ * the serialized expression.
+ */
+ template <class Archive>
+ void
+ serialize(Archive &archive, const unsigned int version);
+
+ //@}
+
+ /**
+ * @name Values
+ */
+ //@{
+
+ /**
+ * Return the value or expression that this class instance represents.
+ */
+ const SE::Expression &
+ get_expression() const;
+
+ /**
+ * Return the primitive SymEngine data type that stores the value or
+ * expression represented by this object.
+ */
+ const SE::Basic &
+ get_value() const;
+
+ /**
+ * Return the pointer to the primitive SymEngine data type that stores
+ * the value or expression represented by this object.
+ */
+ const SE::RCP<const SE::Basic> &
+ get_RCP() const;
+
+ //@}
+
+ /**
+ * @name Math and relational operators with (potentially) symbolic types
+ */
+ //@{
+
+ /**
+ * Assignment operator.
+ *
+ * Sets the data of this object's @p expression equal
+ * to that of the @p rhs object.
+ */
+ Expression &
+ operator=(const Expression &rhs);
+
+ /**
+ * Assignment operator.
+ *
+ * Sets the data of this object's @p expression equal
+ * to that of the @p rhs object.
+ */
+ Expression &
+ operator=(Expression &&rhs);
+
+ /**
+ * Addition assignment.
+ *
+ * The @p rhs @p expression is added in-place to that of
+ * this object's @p expression.
+ */
+ Expression &
+ operator+=(const Expression &rhs);
+
+ /**
+ * Subtraction assignment.
+ *
+ * The @p rhs @p expression is subtracted in-place from that of
+ * this object's @p expression.
+ */
+ Expression &
+ operator-=(const Expression &rhs);
+
+ /**
+ * Multiplication assignment.
+ *
+ * This object's @p expression is multiplied in-place by that of
+ * the @p rhs @p expression.
+ */
+ Expression &
+ operator*=(const Expression &rhs);
+
+ /**
+ * Division assignment.
+ *
+ * This object's @p expression is divided in-place by that of
+ * the @p rhs @p expression.
+ */
+ Expression &
+ operator/=(const Expression &rhs);
+
+ //@}
+
+ /**
+ * @name Math and relational operators with numeric types
+ */
+ //@{
+
+ /**
+ * Assignment operator.
+ *
+ * Set the data of this object's @p expression equal
+ * to the numerical value of the @p rhs.
+ */
+ template <typename NumberType>
+ Expression &
+ operator=(const NumberType &rhs);
+
+ /**
+ * Negation operator.
+ *
+ * Return a the result of pre-multipying this object's @p expression
+ * by <tt>-1</tt>.
+ *
+ * @note This operation is not performed in-place.
+ */
+ Expression
+ operator-() const;
+
+ /**
+ * Addition assignment.
+ *
+ * The numerical value of the @p rhs is added in-place to that of
+ * this object's @p expression.
+ */
+ template <typename NumberType>
+ Expression &
+ operator+=(const NumberType &rhs);
+
+ /**
+ * Subtraction assignment.
+ *
+ * The numerical value of the @p rhs is subtracted in-place from that of
+ * this object's @p expression.
+ */
+ template <typename NumberType>
+ Expression &
+ operator-=(const NumberType &rhs);
+
+ /**
+ * Multiplication assignment.
+ *
+ * This object's @p expression is multiplied in-place by that of
+ * the numerical value of the @p rhs.
+ */
+ template <typename NumberType>
+ Expression &
+ operator*=(const NumberType &rhs);
+
+ /**
+ * Division assignment.
+ *
+ * This object's @p expression is divided in-place by that of
+ * the numerical value of the @p rhs.
+ */
+ template <typename NumberType>
+ Expression &
+ operator/=(const NumberType &rhs);
+
+ //@}
+
+ /**
+ * @name Differentiation
+ */
+ //@{
+
+ /**
+ * Return the derivative of this object's @p expression
+ * with respect to the given @p symbol.
+ */
+ Expression
+ differentiate(const Expression &symbol) const;
+
+ /**
+ * Return the derivative of this object's @p expression
+ * with respect to the given @p symbol.
+ */
+ Expression
+ differentiate(const SE::RCP<const SE::Symbol> &symbol) const;
+
+ /**
+ * Return the derivative of this object's @p expression
+ * with respect to the potential @p symbol.
+ */
+ Expression
+ differentiate(const SE::RCP<const SE::Basic> &symbol) const;
+
+ //@}
+
+ /**
+ * @name Dictionary-based substitution
+ */
+ //@{
+
+ /**
+ * Perform substitution of all symbols found in this object's @p expression
+ * that match a key in the @p substitution_values map.
+ *
+ * @note The replacement value (the entry in the @p substitution_values
+ * that is paired with a key) need not necessarily be numerical, but may
+ * also be another symbolic type.
+ *
+ * @note With dictionary substitution, partial substitution is allowed
+ * (i.e. an incomplete substitution map can be used and the return type
+ * can be symbolic).
+ */
+ Expression
+ substitute(const types::substitution_map &substitution_values) const;
+
+ /**
+ * Perform substitution of all symbols found in this object's @p expression
+ * that match the @p symbol. Each @p symbol will be substituted with
+ * the given @p value.
+ *
+ * @note With dictionary substitution, partial substitution is allowed
+ * (i.e. an incomplete substitution map can be used and the return type
+ * can be symbolic).
+ */
+ Expression
+ substitute(const Expression &symbol, const Expression &value) const;
+
+ /**
+ * Perform substitution of all symbols found in this object's @p expression
+ * that match the @p symbol. Each @p symbol will be substituted with
+ * the given @p value.
+ *
+ * @note With dictionary substitution, partial substitution is allowed
+ * (i.e. an incomplete substitution map can be used and the return type
+ * can be symbolic).
+ */
+ template <typename NumberType>
+ Expression
+ substitute(const Expression &symbol, const NumberType &value) const;
+
+ /**
+ * Full substitution and evaluation. This creates a Expression by
+ * symbol substitution and then immediately computes its numerical value.
+ *
+ * @note All symbols must be resolved by the substitution map in order
+ * for this function to return successfully.
+ */
+ template <typename ReturnType>
+ ReturnType
+ substitute_and_evaluate(
+ const types::substitution_map &substitution_values) const;
+
+ //@}
+
+ /**
+ * @name Conversion operators
+ */
+ //@{
+
+ /**
+ * Conversion operator for real integer or floating point values, and
+ * complex integer or floating point values.
+ *
+ * @note This function is marked explicit so that the conversion must be
+ * performed using a static_cast. In normal use, one would have expected
+ * (Expression)*(double) --> (Expression)
+ * If this function were not marked as explicit, then we could potentially
+ * have
+ * (Expression)*(double) --> (double)
+ * So, to get out a value on needs to do the following:
+ *
+ * <code>
+ * const NumberType val = static_cast<NumberType>(Expression);
+ * </code>
+ *
+ * or, probably less desirably,
+ *
+ * <code>
+ * const NumberType val = NumberType(Expression);
+ * </code>
+ *
+ * @note If the underlying number is a custom type (i.e. encapsulated by a
+ * NumberWrapper), then it is necessary to derive a new class from
+ * Expression and define a specialized conversion operator
+ * that calls an Evaluator that is specialized for this custom number
+ * type. This could be achieved with an overriding conversion function
+ * in the base class, for example:
+ *
+ * <code>
+ * class MyNumber : public Expression
+ * {
+ * ...
+ * template <typename ResultType>
+ * explicit operator ResultType() const
+ * {
+ * if (this->get_value()->get_type_code() == SE::NUMBER_WRAPPER)
+ * {
+ * // Implement custom evaluation function
+ * const ResultType result = ...;
+ * return result;
+ * }
+ * else
+ * // Call base class conversion operator
+ * return Expression::operator ResultType();
+ * }
+ * }
+ * </code>
+ */
+ template <typename ResultType>
+ explicit operator ResultType() const;
+
+ /**
+ * Conversion operator that returns the value or expression that this
+ * class instance represents.
+ */
+ explicit operator const SE::Expression &() const;
+
+ /**
+ * Conversion operator that returns a SymEngine reference counted pointer
+ * to the fundamental type.
+ */
+ operator const SE::RCP<const SE::Basic> &() const;
+
+ //@}
+
+ protected:
+ /**
+ * Return the value or expression that this class instance represents.
+ */
+ SE::Expression &
+ get_expression();
+
+ private:
+ /**
+ * The value or expression that this instance of this class is to
+ * represent.
+ */
+ SE::Expression expression;
+ };
+
+ /**
+ * @name Type traits
+ */
+ //@{
+
+ /**
+ * A struct to indicate whether a given @p NumberType is a supported
+ * symbolically differentiable number or not.
+ * This is a specialization for the deal.II Expression class.
+ *
+ * @author Jean-Paul Pelteret, 2019
+ */
+ template <>
+ struct is_symengine_number<Expression> : std::true_type
+ {};
+
+
+ /**
+ * A struct to indicate whether a given @p NumberType is a supported
+ * SymEngine number or not.
+ * This is a specialization for the deal.II Expression class.
+ *
+ * @author Jean-Paul Pelteret, 2019
+ */
+ template <>
+ struct is_sd_number<Expression> : std::true_type
+ {};
+
+ //@}
+
+ /**
+ * @name Bitwise operators
+ */
+ //@{
+
+ /**
+ * Bitwise left shift operator.
+ *
+ * This is used to output the @p expression to the input @p stream.
+ */
+ std::ostream &
+ operator<<(std::ostream &stream, const Expression &expression);
+
+ /**
+ * Bitwise right shift operator.
+ *
+ * This is used to read in an @p expression from the input @p stream.
+ */
+ std::istream &
+ operator>>(std::istream &stream, Expression &expression);
+
+ //@}
+
+ /**
+ * @name Comparison operators
+ */
+ //@{
+
+ /**
+ * Equality operator.
+ *
+ * Return whether the @p lhs is equal to the @p rhs.
+ */
+ Expression
+ operator==(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * Non-equality operator.
+ *
+ * Return whether the @p lhs is not equal to the @p rhs.
+ */
+ Expression
+ operator!=(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * Less than operator.
+ *
+ * Return whether the @p lhs is less than the @p rhs.
+ */
+ Expression
+ operator<(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * Greater than operator.
+ *
+ * Return whether the @p lhs is greater than the @p rhs.
+ */
+ Expression
+ operator>(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * Less than or equals operator.
+ *
+ * Return whether the @p lhs is less than, or equal to, the @p rhs.
+ */
+ Expression
+ operator<=(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * Greater than or equals operator.
+ *
+ * Return whether the @p lhs is greater than, or equal to, the @p rhs.
+ */
+ Expression
+ operator>=(const Expression &lhs, const Expression &rhs);
+
+ //@}
+
+ /**
+ * @name Logical operators
+ */
+ //@{
+
+ /**
+ * Logical not operator.
+ *
+ * @note This operator can only be applied on boolean and conditional
+ * expressions.
+ */
+ Expression operator!(const Expression &expression);
+
+ /**
+ * Logical and operator.
+ *
+ * @note This operator can only be applied on boolean and conditional
+ * expressions.
+ */
+ Expression operator&(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * Logical inclusive or operator.
+ *
+ * @note This operator can only be applied on boolean and conditional
+ * expressions.
+ */
+ Expression
+ operator|(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * Logical exclusive or (xor) operator.
+ *
+ * @note This operator can only be applied on boolean and conditional
+ * expressions.
+ */
+ Expression
+ operator^(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * And operator.
+ *
+ * @note This operator can only be applied on boolean and conditional
+ * expressions.
+ * This operator is a convenience wrapper for the logical and operator.
+ */
+ Expression
+ operator&&(const Expression &lhs, const Expression &rhs);
+
+ /**
+ * Inclusive or operator.
+ *
+ * @note This operator can only be applied on boolean and conditional
+ * expressions.
+ * This operator is a convenience wrapper for the logical or operator.
+ */
+ Expression
+ operator||(const Expression &lhs, const Expression &rhs);
+
+ //@}
+
+ /**
+ * @name Mathematical operators
+ */
+ //@{
+
+ /**
+ * Addition operator.
+ *
+ * Return the result of adding the @p rhs to the @p lhs.
+ */
+ Expression
+ operator+(Expression lhs, const Expression &rhs);
+
+ /**
+ * Subtraction operator.
+ *
+ * Return the result of subtracting the @p rhs from the @p lhs.
+ */
+ Expression
+ operator-(Expression lhs, const Expression &rhs);
+
+ /**
+ * Multiplication operator.
+ *
+ * Return the result of multiplying the @p lhs by the @p rhs.
+ */
+ Expression operator*(Expression lhs, const Expression &rhs);
+
+ /**
+ * Division operator.
+ *
+ * Return the result of dividing the @p lhs by the @p rhs.
+ */
+ Expression
+ operator/(Expression lhs, const Expression &rhs);
+
+ /**
+ * General addition operator.
+ *
+ * Return the result of adding the @p rhs to the @p lhs.
+ * The @p lhs type may be any supported number type, and the result
+ * is promoted to a Expression. The type conversion makes writing
+ * scalar expressions using Expression more natural.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_constructible<Expression, NumberType>::value>::type>
+ inline Expression
+ operator+(const NumberType &lhs, const Expression &rhs)
+ {
+ return Expression(lhs) + rhs;
+ }
+
+ /**
+ * General addition operator.
+ *
+ * Return the result of adding the @p rhs to the @p lhs.
+ * The @p rhs type may be any supported number type, and the result
+ * is promoted to a Expression. The type conversion makes writing
+ * scalar expressions using Expression more natural.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_constructible<Expression, NumberType>::value>::type>
+ inline Expression
+ operator+(const Expression &lhs, const NumberType &rhs)
+ {
+ return lhs + Expression(rhs);
+ }
+
+ /**
+ * General subtraction operator.
+ *
+ * Return the result of subtracting the @p rhs from the @p lhs.
+ * The @p lhs type may be any supported number type, and the result
+ * is promoted to a Expression. The type conversion makes writing
+ * scalar expressions using Expression more natural.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_constructible<Expression, NumberType>::value>::type>
+ inline Expression
+ operator-(const NumberType &lhs, const Expression &rhs)
+ {
+ return Expression(lhs) - rhs;
+ }
+
+ /**
+ * General subtraction operator.
+ *
+ * Return the result of subtracting the @p rhs from the @p lhs.
+ * The @p rhs type may be any supported number type, and the result
+ * is promoted to a Expression. The type conversion makes writing
+ * scalar expressions using Expression more natural.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_constructible<Expression, NumberType>::value>::type>
+ inline Expression
+ operator-(const Expression &lhs, const NumberType &rhs)
+ {
+ return lhs - Expression(rhs);
+ }
+
+ /**
+ * General multiplication operator.
+ *
+ * Return the result of multiplying the @p lhs by the @p rhs.
+ * The @p lhs type may be any supported number type, and the result
+ * is promoted to a Expression. The type conversion makes writing
+ * scalar expressions using Expression more natural.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_constructible<Expression, NumberType>::value>::type>
+ inline Expression operator*(const NumberType &lhs, const Expression &rhs)
+ {
+ return Expression(lhs) * rhs;
+ }
+
+ /**
+ * General multiplication operator.
+ *
+ * Return the result of multiplying the @p lhs by the @p rhs.
+ * The @p rhs type may be any supported number type, and the result
+ * is promoted to a Expression. The type conversion makes writing
+ * scalar expressions using Expression more natural.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_constructible<Expression, NumberType>::value>::type>
+ inline Expression operator*(const Expression &lhs, const NumberType &rhs)
+ {
+ return lhs * Expression(rhs);
+ }
+
+ /**
+ * General division operator.
+ *
+ * Return the result of dividing the @p lhs by the @p rhs.
+ * The @p lhs type may be any supported number type, and the result
+ * is promoted to a Expression. The type conversion makes writing
+ * scalar expressions using Expression more natural.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_constructible<Expression, NumberType>::value>::type>
+ inline Expression
+ operator/(const NumberType &lhs, const Expression &rhs)
+ {
+ return Expression(lhs) / rhs;
+ }
+
+ /**
+ * General division operator.
+ *
+ * Return the result of dividing the @p lhs by the @p rhs.
+ * The @p lhs type may be any supported number type, and the result
+ * is promoted to a Expression. The type conversion makes writing
+ * scalar expressions using Expression more natural.
+ */
+ template <typename NumberType,
+ typename = typename std::enable_if<
+ std::is_constructible<Expression, NumberType>::value>::type>
+ inline Expression
+ operator/(const Expression &lhs, const NumberType &rhs)
+ {
+ return lhs / Expression(rhs);
+ }
+
+ //@}
+
+ } // namespace SD
+} // namespace Differentiation
+
+
+/* ---------------------- inline and template functions -------------------- */
+
+
+# ifndef DOXYGEN
+
+
+namespace Differentiation
+{
+ namespace SD
+ {
+ template <typename NumberType, typename>
+ Expression::Expression(const NumberType &value)
+ : expression(value)
+ {}
+
+
+ template <typename NumberType, typename>
+ Expression::Expression(const std::complex<NumberType> &value)
+ : expression(value)
+ {}
+
+
+ template <typename NumberType, typename>
+ Expression::Expression(const NumberType &numerator,
+ const NumberType &denominator)
+ : expression(SE::Rational::from_two_ints(*SE::integer(numerator),
+ *SE::integer(denominator)))
+ {}
+
+
+ template <class Archive>
+ void
+ Expression::serialize(Archive &ar, const unsigned int /*version*/)
+ {
+ // This is a bit tricky... The SymEngine expression class is not
+ // serializable, and we're not sure if we're writing out or
+ // reading in here. So what we'll do is try to synchronise the
+ // current data with the stream, predicting if we're likely trying
+ // to write out data. We do this by checking if we're currently working
+ // with an object that has a trivial state. If not, then we're likely
+ // to output this object's contents.
+ const Expression default_constructed;
+ const bool likely_writing_out =
+ (get_RCP()->__eq__(*default_constructed.get_RCP()) == false);
+
+ std::stringstream sstream;
+ std::string expr = default_constructed.get_RCP()->__str__();
+
+ // Output
+ if (likely_writing_out)
+ {
+ sstream << *this;
+ expr = sstream.str();
+ }
+
+ // Serialise
+ ar &expr;
+
+ // Input
+ if (!likely_writing_out)
+ {
+ sstream.clear();
+ sstream << expr;
+ sstream >> *this;
+
+ parse(expr);
+ }
+ }
+
+
+ template <typename NumberType>
+ Expression
+ Expression::substitute(const Expression &symbol,
+ const NumberType &value) const
+ {
+ Assert(SE::is_a<SE::Symbol>(symbol.get_value()),
+ ExcMessage(
+ "Substitution with a number that does not represent a symbol."));
+
+ types::substitution_map sub_vals;
+ sub_vals[symbol] = Expression(value);
+ return substitute(sub_vals);
+ }
+
+
+ template <typename ReturnType>
+ ReturnType
+ Expression::substitute_and_evaluate(
+ const types::substitution_map &substitution_values) const
+ {
+ return static_cast<ReturnType>(substitute(substitution_values));
+ }
+
+
+ template <typename NumberType>
+ Expression &
+ Expression::operator=(const NumberType &rhs)
+ {
+ *this = Expression(rhs);
+ return *this;
+ }
+
+
+ template <typename NumberType>
+ Expression &
+ Expression::operator+=(const NumberType &rhs)
+ {
+ *this = Expression(SE::add(get_RCP(), Expression(rhs).get_RCP()));
+ return *this;
+ }
+
+
+ template <typename NumberType>
+ Expression &
+ Expression::operator-=(const NumberType &rhs)
+ {
+ *this = Expression(SE::sub(get_RCP(), Expression(rhs).get_RCP()));
+ return *this;
+ }
+
+
+ template <typename NumberType>
+ Expression &
+ Expression::operator*=(const NumberType &rhs)
+ {
+ *this = Expression(SE::mul(get_RCP(), Expression(rhs).get_RCP()));
+ return *this;
+ }
+
+
+ template <typename NumberType>
+ Expression &
+ Expression::operator/=(const NumberType &rhs)
+ {
+ *this = Expression(SE::div(get_RCP(), Expression(rhs).get_RCP()));
+ return *this;
+ }
+
+
+ template <typename ResultType>
+ Expression::operator ResultType() const
+ {
+ return static_cast<ResultType>(get_expression());
+ }
+
+ } // namespace SD
+} // namespace Differentiation
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+# endif // DOXYGEN
+
+#endif // DEAL_II_WITH_SYMENGINE
+
+#endif // dealii_differentiation_sd_symengine_number_types_h
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SYMENGINE
+
+# include <deal.II/base/signaling_nan.h>
+
+# include <deal.II/differentiation/sd/symengine_number_types.h>
+# include <deal.II/differentiation/sd/symengine_types.h>
+# include <deal.II/differentiation/sd/symengine_utilities.h>
+
+# include <symengine/complex_double.h>
+# include <symengine/logic.h>
+# include <symengine/number.h>
+# include <symengine/parser.h>
+# include <symengine/real_double.h>
+# include <symengine/symbol.h>
+# include <symengine/symengine_exception.h>
+
+# include <string>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Differentiation
+{
+ namespace SD
+ {
+ namespace SE = ::SymEngine;
+
+
+ /* ---------------------------- Constructors -------------------------- */
+
+
+ Expression::Expression()
+ : expression()
+ {}
+
+
+ Expression::Expression(const bool &value)
+ : expression(SE::boolean(value))
+ {}
+
+
+ Expression::Expression(const SE::integer_class &value)
+ : expression(value)
+ {}
+
+
+ Expression::Expression(const SE::rational_class &value)
+ : expression(value)
+ {}
+
+
+ Expression::Expression(const Expression &condition,
+ const Expression &expression_if_true,
+ const Expression &expression_if_false)
+ {
+ Assert(SE::is_a_Boolean(condition.get_value()),
+ ExcMessage(
+ "The conditional expression must return a boolean type."));
+
+ const SE::RCP<const SE::Boolean> condition_rcp =
+ SE::rcp_static_cast<const SE::Boolean>(condition.get_RCP());
+ expression =
+ SE::piecewise({{expression_if_true.get_RCP(), condition_rcp},
+ {expression_if_false.get_RCP(), SE::boolTrue}});
+ }
+
+
+ Expression::Expression(const std::vector<std::pair<Expression, Expression>>
+ & condition_expression,
+ const Expression &expression_otherwise)
+ {
+ SE::PiecewiseVec piecewise_function;
+ piecewise_function.reserve(condition_expression.size() + 1);
+
+ // Add tested conditional entries
+ for (const auto &entry : condition_expression)
+ {
+ Assert(SE::is_a_Boolean(entry.first.get_value()),
+ ExcMessage(
+ "The conditional expression must return a boolean type."));
+ piecewise_function.push_back(
+ {entry.second.get_RCP(),
+ SE::rcp_static_cast<const SE::Boolean>(entry.first.get_RCP())});
+ }
+
+ // Add default value
+ piecewise_function.push_back(
+ {expression_otherwise.get_RCP(), SE::boolTrue});
+
+ // Initialize
+ expression = SE::piecewise(std::move(piecewise_function));
+ }
+
+
+ Expression::Expression(const std::vector<std::pair<Expression, Expression>>
+ &condition_expression)
+ {
+ // Use the other constructor with a fatal termination point
+ // ensuring that an error is thrown if none of the conditions
+ // are met.
+ *this = Expression(condition_expression,
+ Expression(numbers::signaling_nan<double>()));
+ }
+
+
+ Expression::Expression(const char *symbol)
+ : expression(SE::symbol(symbol))
+ {}
+
+
+ Expression::Expression(const std::string &str,
+ const bool parse_as_expression)
+ {
+ try
+ {
+ expression = (parse_as_expression ?
+ SE::parse(str) // The string is a symbolic "name"
+ :
+ SE::rcp_static_cast<const SE::Basic>(SE::symbol(
+ str))); // The string is a symbolic "expression"
+ }
+ catch (...)
+ {
+ AssertThrow(false, ExcSymEngineParserError(str));
+ }
+ }
+
+
+ Expression::Expression(const std::string & symbol_func,
+ const types::symbol_vector &arguments)
+ : expression(SE::function_symbol(
+ symbol_func,
+ Utilities::convert_expression_vector_to_basic_vector(arguments)))
+ {}
+
+
+ Expression::Expression(const SE::Expression &rhs)
+ : expression(rhs)
+ {}
+
+
+ Expression::Expression(const SE::RCP<const SE::Basic> &rhs)
+ : expression(rhs)
+ {}
+
+
+ Expression::Expression(SE::RCP<const SE::Basic> &&rhs)
+ : expression(rhs)
+ {}
+
+
+ /* ------------------------------ Utilities ---------------------------- */
+
+
+ Expression &
+ Expression::parse(const std::string &expression)
+ {
+ *this = Expression(expression, true /*parse_as_expression*/);
+ return *this;
+ }
+
+
+ std::ostream &
+ Expression::print(std::ostream &os) const
+ {
+ os << *this;
+ return os;
+ }
+
+
+ void
+ Expression::save(std::ostream &os) const
+ {
+ // We write each expression on a new line.
+ // Note: SymEngine outputs a non-terminating string
+ os << *this;
+ os << std::endl;
+ }
+
+
+ void
+ Expression::load(std::istream &is)
+ {
+ // Need to make sure that we read the entire line in,
+ // and then subsequently parse it.
+ std::string expr;
+ std::getline(is, expr);
+ Assert(!is.bad(), ExcIO());
+ parse(expr);
+ }
+
+
+ /* ------------------------------- Values ----------------------------- */
+
+
+ const SE::Expression &
+ Expression::get_expression() const
+ {
+ return expression;
+ }
+
+
+ SE::Expression &
+ Expression::get_expression()
+ {
+ return expression;
+ }
+
+
+ const SE::Basic &
+ Expression::get_value() const
+ {
+ return *get_RCP();
+ }
+
+
+ const SE::RCP<const SE::Basic> &
+ Expression::get_RCP() const
+ {
+ return get_expression().get_basic();
+ }
+
+
+ /* --------------------------- Differentiation ------------------------- */
+
+
+ Expression
+ Expression::differentiate(const SE::RCP<const SE::Symbol> &symbol) const
+ {
+ return Expression(SE::diff(get_RCP(), symbol));
+ }
+
+
+ Expression
+ Expression::differentiate(const SE::RCP<const SE::Basic> &symbol) const
+ {
+ // Potential symbol
+ return Expression(SE::sdiff(get_RCP(), symbol));
+ }
+
+
+ Expression
+ Expression::differentiate(const Expression &symbol) const
+ {
+ return differentiate(symbol.get_RCP());
+ }
+
+
+ /* ------------- Conversion operators ------------------------- */
+
+
+ Expression::operator const SE::Expression &() const
+ {
+ return get_expression();
+ }
+
+
+ Expression::operator const SE::RCP<const SE::Basic> &() const
+ {
+ return get_expression().get_basic();
+ }
+
+
+ /* ------------- Dictionary-based substitution ------------------------- */
+
+
+ Expression
+ Expression::substitute(
+ const types::substitution_map &substitution_values) const
+ {
+ return Expression(get_expression().subs(
+ Utilities::convert_expression_map_to_basic_map(substitution_values)));
+ }
+
+
+ Expression
+ Expression::substitute(const Expression &symbol,
+ const Expression &value) const
+ {
+ Assert(SE::is_a<SE::Symbol>(symbol.get_value()),
+ ExcMessage(
+ "Substitution with a number that does not represent a symbol."));
+
+ types::substitution_map sub_vals;
+ sub_vals[symbol] = value;
+ return substitute(sub_vals);
+ }
+
+
+ /* -------------------- Math and relational operators ------------------ */
+
+
+ Expression &
+ Expression::operator=(const Expression &rhs)
+ {
+ if (this != &rhs)
+ this->expression = rhs.get_expression();
+
+ return *this;
+ }
+
+
+ Expression &
+ Expression::operator=(Expression &&rhs)
+ {
+ if (this != &rhs)
+ this->expression = std::move(rhs.expression);
+
+ return *this;
+ }
+
+
+ Expression
+ Expression::operator-() const
+ {
+ return Expression(-get_expression());
+ }
+
+
+ Expression &
+ Expression::operator+=(const Expression &rhs)
+ {
+ this->expression += rhs.get_expression();
+ return *this;
+ }
+
+
+ Expression &
+ Expression::operator-=(const Expression &rhs)
+ {
+ this->expression -= rhs.get_expression();
+ return *this;
+ }
+
+
+ Expression &
+ Expression::operator*=(const Expression &rhs)
+ {
+ this->expression *= rhs.get_expression();
+ return *this;
+ }
+
+
+ Expression &
+ Expression::operator/=(const Expression &rhs)
+ {
+ this->expression /= rhs.get_expression();
+ return *this;
+ }
+
+
+ std::ostream &
+ operator<<(std::ostream &stream, const Expression &expr)
+ {
+ stream << expr.get_expression();
+ return stream;
+ }
+
+
+ std::istream &
+ operator>>(std::istream &stream, Expression &expr)
+ {
+ std::string str;
+ stream >> str;
+ expr.parse(str);
+ return stream;
+ }
+
+
+ Expression
+ operator==(const Expression &lhs, const Expression &rhs)
+ {
+ return Expression(SE::Eq(lhs.get_RCP(), rhs.get_RCP()));
+ }
+
+
+ Expression
+ operator!=(const Expression &lhs, const Expression &rhs)
+ {
+ return Expression(SE::Ne(lhs.get_RCP(), rhs.get_RCP()));
+ }
+
+
+ Expression
+ operator<(const Expression &lhs, const Expression &rhs)
+ {
+ return Expression(SE::Lt(lhs.get_RCP(), rhs.get_RCP()));
+ }
+
+
+ Expression
+ operator>(const Expression &lhs, const Expression &rhs)
+ {
+ return Expression(SE::Gt(lhs.get_RCP(), rhs.get_RCP()));
+ }
+
+
+ Expression
+ operator<=(const Expression &lhs, const Expression &rhs)
+ {
+ return Expression(SE::Le(lhs.get_RCP(), rhs.get_RCP()));
+ }
+
+
+ Expression
+ operator>=(const Expression &lhs, const Expression &rhs)
+ {
+ return Expression(SE::Ge(lhs.get_RCP(), rhs.get_RCP()));
+ }
+
+
+ Expression operator!(const Expression &expression)
+ {
+ Assert(SE::is_a_Boolean(expression.get_value()),
+ ExcMessage("The expression must return a boolean type."));
+
+ const SE::RCP<const SE::Boolean> expression_rcp =
+ SE::rcp_static_cast<const SE::Boolean>(expression.get_RCP());
+
+ return Expression(SE::logical_not(expression_rcp));
+ }
+
+
+ Expression operator&(const Expression &lhs, const Expression &rhs)
+ {
+ Assert(SE::is_a_Boolean(lhs.get_value()),
+ ExcMessage("The lhs expression must return a boolean type."));
+ Assert(SE::is_a_Boolean(rhs.get_value()),
+ ExcMessage("The rhs expression must return a boolean type."));
+
+ const SE::RCP<const SE::Boolean> lhs_rcp =
+ SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
+ const SE::RCP<const SE::Boolean> rhs_rcp =
+ SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
+
+ return Expression(SE::logical_and({lhs_rcp, rhs_rcp}));
+ }
+
+
+ Expression
+ operator|(const Expression &lhs, const Expression &rhs)
+ {
+ Assert(SE::is_a_Boolean(lhs.get_value()),
+ ExcMessage("The lhs expression must return a boolean type."));
+ Assert(SE::is_a_Boolean(rhs.get_value()),
+ ExcMessage("The rhs expression must return a boolean type."));
+
+ const SE::RCP<const SE::Boolean> lhs_rcp =
+ SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
+ const SE::RCP<const SE::Boolean> rhs_rcp =
+ SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
+
+ return Expression(SE::logical_or({lhs_rcp, rhs_rcp}));
+ }
+
+
+ Expression
+ operator^(const Expression &lhs, const Expression &rhs)
+ {
+ Assert(SE::is_a_Boolean(lhs.get_value()),
+ ExcMessage("The lhs expression must return a boolean type."));
+ Assert(SE::is_a_Boolean(rhs.get_value()),
+ ExcMessage("The rhs expression must return a boolean type."));
+
+ const SE::RCP<const SE::Boolean> lhs_rcp =
+ SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
+ const SE::RCP<const SE::Boolean> rhs_rcp =
+ SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
+
+ return Expression(SE::logical_xor({lhs_rcp, rhs_rcp}));
+ }
+
+
+ Expression
+ operator&&(const Expression &lhs, const Expression &rhs)
+ {
+ return lhs & rhs;
+ }
+
+
+ Expression
+ operator||(const Expression &lhs, const Expression &rhs)
+ {
+ return lhs | rhs;
+ }
+
+
+ Expression
+ operator+(Expression lhs, const Expression &rhs)
+ {
+ lhs += rhs;
+ return lhs;
+ }
+
+
+ Expression
+ operator-(Expression lhs, const Expression &rhs)
+ {
+ lhs -= rhs;
+ return lhs;
+ }
+
+
+ Expression operator*(Expression lhs, const Expression &rhs)
+ {
+ lhs *= rhs;
+ return lhs;
+ }
+
+
+ Expression
+ operator/(Expression lhs, const Expression &rhs)
+ {
+ lhs /= rhs;
+ return lhs;
+ }
+
+
+ } // namespace SD
+} // namespace Differentiation
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SYMENGINE
projects / dealii.git / commitdiff