]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add basis for SymEngine (SD) wrappers
authorJean-Paul Pelteret <jppelteret@gmail.com>
Sat, 30 Mar 2019 10:19:55 +0000 (11:19 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Tue, 2 Apr 2019 16:38:38 +0000 (18:38 +0200)
- Types
- Expression class
- Utilities

doc/news/changes/major/20190329Jean-PaulPelteret [new file with mode: 0644]
include/deal.II/differentiation/sd/symengine_number_types.h [new file with mode: 0644]
include/deal.II/differentiation/sd/symengine_types.h [new file with mode: 0644]
include/deal.II/differentiation/sd/symengine_utilities.h [new file with mode: 0644]
source/differentiation/CMakeLists.txt
source/differentiation/sd/CMakeLists.txt [new file with mode: 0644]
source/differentiation/sd/symengine_number_types.cc [new file with mode: 0644]
source/differentiation/sd/symengine_types.cc [new file with mode: 0644]
source/differentiation/sd/symengine_utilities.cc [new file with mode: 0644]

diff --git a/doc/news/changes/major/20190329Jean-PaulPelteret b/doc/news/changes/major/20190329Jean-PaulPelteret
new file mode 100644 (file)
index 0000000..d5f571d
--- /dev/null
@@ -0,0 +1,16 @@
+New: A new class Differentiation::SD::Expression has been added. which allows
+the creation and manipulation of scalar symbolic functions using the SymEngine
+library. It will later be used as the basis of symbolic tensor calculus, amongst
+other things. It incorporates numerous features, which currently include:
+<ul>
+  <li>expression parsing,</li>
+  <li>comparison operations,</li>
+  <li>logical operations,</li>
+  <li>basic math operations,</li>
+  <li>conditional expression construction,</li>
+  <li>differentiation,</li>
+  <li>substitution (partial and complete), and</li>
+  <li>serialization.</li>
+</ul>
+<br>
+(Jean-Paul Pelteret, 2019/03/29)
diff --git a/include/deal.II/differentiation/sd/symengine_number_types.h b/include/deal.II/differentiation/sd/symengine_number_types.h
new file mode 100644 (file)
index 0000000..b3a487d
--- /dev/null
@@ -0,0 +1,1302 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/include/deal.II/differentiation/sd/symengine_types.h b/include/deal.II/differentiation/sd/symengine_types.h
new file mode 100644 (file)
index 0000000..1f493c0
--- /dev/null
@@ -0,0 +1,81 @@
+// ---------------------------------------------------------------------
+//
+// 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_types_h
+#define dealii_differentiation_sd_symengine_types_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SYMENGINE
+
+#  include <map>
+#  include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Differentiation
+{
+  namespace SD
+  {
+    // Forward declarations
+    class Expression;
+
+
+    namespace types
+    {
+      namespace internal
+      {
+        /**
+         * A comparator for Expressions used as keys in maps.
+         */
+        struct ExpressionKeyLess
+        {
+          bool
+          operator()(const SD::Expression &lhs,
+                     const SD::Expression &rhs) const;
+        };
+      } // namespace internal
+
+      /**
+       * Type definition for a value substitution map.
+       *
+       * This serves the same purpose as a
+       * <code>SymEngine::map_basic_basic</code>, which is equivalent to a
+       * <code>std::map<SE::RCP<const SE::Basic>,  SE::RCP<const
+       * SE::Basic>></code>.
+       */
+      using substitution_map =
+        std::map<SD::Expression, SD::Expression, internal::ExpressionKeyLess>;
+
+      /**
+       * Type definition for a vector of symbols.
+       *
+       * This serves the same purpose as a <code>SymEngine::vec_basic</code>,
+       * which is equivalent to a <code>std::vector<SE::RCP<const
+       * SE::Basic>></code>
+       */
+      using symbol_vector = std::vector<SD::Expression>;
+
+    } // namespace types
+
+  } // namespace SD
+} // namespace Differentiation
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SYMENGINE
+
+#endif // dealii_differentiation_sd_symengine_types_h
diff --git a/include/deal.II/differentiation/sd/symengine_utilities.h b/include/deal.II/differentiation/sd/symengine_utilities.h
new file mode 100644 (file)
index 0000000..e9f64c5
--- /dev/null
@@ -0,0 +1,62 @@
+// ---------------------------------------------------------------------
+//
+// 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_utilities_h
+#define dealii_differentiation_sd_symengine_utilities_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SYMENGINE
+
+#  include <deal.II/differentiation/sd/symengine_types.h>
+
+#  include <symengine/basic.h>
+#  include <symengine/dict.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Differentiation
+{
+  namespace SD
+  {
+    namespace SE = ::SymEngine;
+
+    namespace Utilities
+    {
+      /**
+       * Convert a map of Expressions to its SymEngine counterpart.
+       */
+      SE::map_basic_basic
+      convert_expression_map_to_basic_map(
+        const SD::types::substitution_map &substitution_map);
+
+      /**
+       * Convert a vector of Expressions to its SymEngine counterpart.
+       */
+      SE::vec_basic
+      convert_expression_vector_to_basic_vector(
+        const SD::types::symbol_vector &symbol_vector);
+
+    } // namespace Utilities
+
+  } // namespace SD
+} // namespace Differentiation
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SYMENGINE
+
+#endif // dealii_differentiation_sd_symengine_utilities_h
index 28015cc9d5a74706e727b947806fe1290b03c3fa..3d27c72a33e4b47559a3c7bf32abcd34d3e01776 100644 (file)
@@ -1,6 +1,6 @@
 ## ---------------------------------------------------------------------
 ##
-## Copyright (C) 2017 by the deal.II authors
+## Copyright (C) 2017-2019 by the deal.II authors
 ##
 ## This file is part of the deal.II library.
 ##
@@ -14,3 +14,4 @@
 ## ---------------------------------------------------------------------
 
 ADD_SUBDIRECTORY(ad)
+ADD_SUBDIRECTORY(sd)
diff --git a/source/differentiation/sd/CMakeLists.txt b/source/differentiation/sd/CMakeLists.txt
new file mode 100644 (file)
index 0000000..f673058
--- /dev/null
@@ -0,0 +1,32 @@
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2017 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_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
+
+SET(_src
+  symengine_number_types.cc
+  symengine_types.cc
+  symengine_utilities.cc
+  )
+
+SET(_inst
+  )
+
+FILE(GLOB _header
+  ${CMAKE_SOURCE_DIR}/include/deal.II/differentiation/sd/*.h
+  )
+
+DEAL_II_ADD_LIBRARY(obj_differentiation_sd OBJECT ${_src} ${_header} ${_inst})
+EXPAND_INSTANTIATIONS(obj_differentiation_sd "${_inst}")
diff --git a/source/differentiation/sd/symengine_number_types.cc b/source/differentiation/sd/symengine_number_types.cc
new file mode 100644 (file)
index 0000000..fce4947
--- /dev/null
@@ -0,0 +1,541 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/source/differentiation/sd/symengine_types.cc b/source/differentiation/sd/symengine_types.cc
new file mode 100644 (file)
index 0000000..4c628ae
--- /dev/null
@@ -0,0 +1,55 @@
+// ---------------------------------------------------------------------
+//
+// 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/differentiation/sd/symengine_number_types.h>
+#  include <deal.II/differentiation/sd/symengine_types.h>
+
+#  include <symengine/basic.h>
+#  include <symengine/dict.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Differentiation
+{
+  namespace SD
+  {
+    namespace SE = ::SymEngine;
+
+    namespace types
+    {
+      namespace internal
+      {
+        bool
+        ExpressionKeyLess::operator()(const SD::Expression &lhs,
+                                      const SD::Expression &rhs) const
+        {
+          return SE::RCPBasicKeyLess()(lhs.get_RCP(), rhs.get_RCP());
+        }
+      } // namespace internal
+
+    } // namespace types
+
+  } // namespace SD
+} // namespace Differentiation
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SYMENGINE
diff --git a/source/differentiation/sd/symengine_utilities.cc b/source/differentiation/sd/symengine_utilities.cc
new file mode 100644 (file)
index 0000000..9697316
--- /dev/null
@@ -0,0 +1,63 @@
+// ---------------------------------------------------------------------
+//
+// 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/differentiation/sd/symengine_number_types.h>
+#  include <deal.II/differentiation/sd/symengine_utilities.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Differentiation
+{
+  namespace SD
+  {
+    namespace SE = ::SymEngine;
+
+    namespace Utilities
+    {
+      SE::map_basic_basic
+      convert_expression_map_to_basic_map(
+        const SD::types::substitution_map &substitution_map)
+      {
+        SE::map_basic_basic sub_map;
+        for (const auto &entry : substitution_map)
+          sub_map[entry.first.get_RCP()] = entry.second.get_RCP();
+        return sub_map;
+      }
+
+
+      SE::vec_basic
+      convert_expression_vector_to_basic_vector(
+        const SD::types::symbol_vector &symbol_vector)
+      {
+        SE::vec_basic symb_vec;
+        symb_vec.reserve(symbol_vector.size());
+        for (const auto &entry : symbol_vector)
+          symb_vec.push_back(entry.get_RCP());
+        return symb_vec;
+      }
+
+    } // namespace Utilities
+
+  } // namespace SD
+} // namespace Differentiation
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SYMENGINE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.