]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Symbolic function. 9171/head
authorLuca Heltai <luca.heltai@sissa.it>
Sat, 14 Dec 2019 20:03:03 +0000 (21:03 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 28 Jan 2020 09:19:49 +0000 (10:19 +0100)
15 files changed:
doc/news/changes/major/20191214LucaHeltai [new file with mode: 0644]
include/deal.II/base/function_parser.h
include/deal.II/base/symbolic_function.h [new file with mode: 0644]
include/deal.II/base/symbolic_function.templates.h [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/symbolic_function.cc [new file with mode: 0644]
source/base/symbolic_function.inst.in [new file with mode: 0644]
tests/base/symbolic_function_01.cc [new file with mode: 0644]
tests/base/symbolic_function_01.with_symengine=true.output [new file with mode: 0644]
tests/base/symbolic_function_02.cc [new file with mode: 0644]
tests/base/symbolic_function_02.with_symengine=true.output [new file with mode: 0644]
tests/base/symbolic_function_03.cc [new file with mode: 0644]
tests/base/symbolic_function_03.with_symengine=true.output [new file with mode: 0644]
tests/base/symbolic_function_04.cc [new file with mode: 0644]
tests/base/symbolic_function_04.with_symengine=true.output [new file with mode: 0644]

diff --git a/doc/news/changes/major/20191214LucaHeltai b/doc/news/changes/major/20191214LucaHeltai
new file mode 100644 (file)
index 0000000..5f835f3
--- /dev/null
@@ -0,0 +1,6 @@
+New: The SymbolicFunction<dim> class allows one to leaverage the SymEngine library to generate
+dealii::Function objects where the gradients, Laplacians, and Hessians are computed symbolically
+providing also the possibility to extract the time derivative of the SymbolicFunction<dim> object
+as another SymbolicFunction<dim> object.
+<br>
+(Luca Heltai, 2019/12/14)
index 87a6a5b1d4a3a70a104b0976425ade6fec1a5fc7..1dc24a9cfeba8fdc3db50e6462a12aef31fab5b0 100644 (file)
@@ -211,6 +211,11 @@ class Vector;
  *                        constants);
  * @endcode
  *
+ * @note The difference between this class and the SymbolicFunction class is
+ * that the SymbolicFunction class allows to compute first and second order
+ * derivatives (in a symbolic way), while this class computes first order
+ * derivatives only, using finite differences. For complicated expressions,
+ * this class is generally faster than SymbolicFunction.
  *
  * @ingroup functions
  * @author Luca Heltai, Timo Heister 2005, 2014
diff --git a/include/deal.II/base/symbolic_function.h b/include/deal.II/base/symbolic_function.h
new file mode 100644 (file)
index 0000000..c188d68
--- /dev/null
@@ -0,0 +1,483 @@
+// ---------------------------------------------------------------------
+//
+// 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_symbolic_function_h
+#define dealii_symbolic_function_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/sd.h>
+
+#include <functional>
+#include <iostream>
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+// Forward declarations
+#ifndef DOXYGEN
+template <typename number>
+class Vector;
+namespace Functions
+{
+  template <int dim, typename RangeNumberType>
+  class SymbolicFunction;
+}
+#endif
+
+namespace Functions
+{
+#ifdef DEAL_II_WITH_SYMENGINE
+  /**
+   * A Function class that leverages symbolic differentiation to compute
+   * gradients, Laplacians, Hessians, and time derivatives.
+   *
+   * This class can be used to define functions using methods provided by the
+   * Differentiation::SD namespace. In particular, one can define a symbolic
+   * evaluation point (the argument of the function), as well as a symbolic
+   * expression.
+   *
+   * The symbolic gradients and the symbolic Hessians are computed at
+   * construction time, and when a substitution in the symbolic functions is
+   * requested by the user using the method update_user_substitution_map().
+   *
+   * Whenever one of the evaluation methods is called, a substitution is
+   * attempted with the coordinate symbols argument replaced by the evaluation
+   * point and the symbolic time replaced by the current time, as returned by
+   * the get_time() method. The user has to make sure that at evaluation time
+   * argument substitution provides a fully evaluated expression (i.e., no other
+   * symbols are contained in the function expression, except numerical values),
+   * or an exception will be thrown. Additional symbols can be partially
+   * evaluated or substituted by storing them in a user supplied substitution
+   * maps, that can be updated by calling update_user_substitution_map() or the
+   * set_additional_function_arguments() methods.
+   *
+   * The simplest use case of this class is given in the following example:
+   * @code
+   * SymbolicFunction<2> fun("x^2+y; t*x*y");
+   * fun.set_time(3.0);
+   * Point<2> p(1.0, 2.0);
+   *
+   * auto a = fun.value(p, / * component * / 0); // a = 3.0
+   * auto b = fun.value(p, / * component * / 1); // b = 6.0
+   *
+   * auto df_dt = fun.time_derivative();
+   *
+   * auto c = df_dt.value(p, / * component * / 0); // c = 0.0
+   * auto d = df_dt.value(p, / * component * / 1); // d = 2.0
+   * @endcode
+   * where a Function with two components is defined using a string containing
+   * their expressions separated by semicolons.
+   *
+   * A more involved example, that explicitly uses
+   * Differentiation::SD::Expression objects, is given by
+   * @code
+   * using namespace Differentiation::SD;
+   * // Create a position Tensor<1,2,Differentiation::SD::Expression>
+   * // with symbols "x" and "y", and the symbol "t"
+   * const auto x = SymbolicFunction<2>::get_default_coordinate_symbols();
+   * const auto t = make_symbol("t");
+   *
+   * // Use directly x[0] (the symbol "x"), x[1] (the symbol "y"), and t
+   * // (the symbol "t").
+   * Expression f = std::sin(x[0])*std::cos(x[1])*std::sin(t);
+   * // Alternatively, you can achieve the same result parsing a string:
+   * // Expression f("sin(x)*cos(y)*sin(t)", true);
+   * SymbolicFunction<2> function({f}, x);
+   *
+   * // Evaluate the function, its gradient, and its Laplacian
+   * Point<2> p(1.0, 2.0);
+   * auto fp = function.value(p);
+   * auto gradfp = function.gradient(p);
+   * auto lapfp = function.laplacian(p);
+   *
+   * // Evaluate the time derivative of the function, its gradient, and its
+   * // Laplacian
+   * auto time_derivative = function.time_derivative();
+   * auto dt_fp = time_derivative.value(p);
+   * auto dt_gradfp = time_derivative.gradient(p);
+   * auto dt_lapfp = time_derivative.laplacian(p);
+   * @endcode
+   *
+   * Partial substitution is possible (i.e., you can define the function using
+   * additional symbols). However, as soon as you evaluate the function, you
+   * have to make sure that all extraneous symbols (i.e., those not referring
+   * to the spacial @p coordinate_symbols or to the @p time_symbol variable)
+   * have been substituted with numerical values, or expressions of the spatial
+   * or temporal argument, by calling the update_user_substitution_map() or the
+   * set_additional_function_arguments() methods.
+   *
+   * If your function requires additional arguments to be evaluated, you can
+   * specify them by calling the set_additional_function_arguments() method.
+   *
+   * If you call update_user_substitution_map() and
+   * set_additional_function_arguments() with the same argument, the effect on
+   * the function evaluation will be the same, however, the internal behaviour
+   * and function derivatives will be different. The method
+   * update_user_substitution_map() performs the substitution once (the first
+   * time it is required), and then stores internally a copy of the resulting
+   * expression, together with its derivatives (if required). These are then
+   * used in all subsequent evaluations. Calling
+   * set_additional_function_arguments() will evaluate the passed
+   * subsitution map on the fly during evaluation time, *after* all
+   * derivatives have been computed.
+   *
+   * @note The difference between this class and the FunctionParser class is
+   * that this class allows to compute first and second order derivatives (in a
+   * symbolic way), while the FunctionParser class computes first order
+   * derivatives only, using finite differences. For complicated expressions,
+   * this class may be slower than the FunctionParser class.
+   *
+   * @ingroup functions
+   * @author Luca Heltai 2019
+   */
+  template <int dim, typename RangeNumberType = double>
+  class SymbolicFunction : public Function<dim, RangeNumberType>
+  {
+  public:
+    /**
+     * Constructor.
+     *
+     * The resulting Function object will have as many components as there
+     * are entries in the vector of symbolic expressions @p function.
+     *
+     * The vector @p function should contain a list of symbolic expression
+     * involving the coordinate symbols argument @p coordinate_symbols and
+     * possibly the symbolic time argument @p time_symbol. It is possible to
+     * define it in terms of other symbols, as long as the optional parameter
+     * @p user_substitution_map replaces all symbols except
+     * @p coordinate_symbols and @p time_symbol.
+     * This is useful if, for example, you want to express formulas in terms of
+     * material parameters that you want to name symbolically, rather than
+     * through their numeric values when defining the formula, or when you want
+     * to express your formula in terms of polar coordinates rather than
+     * cartesian ones, and you want the symbolic engine to compute the
+     * derivatives for you.
+     * You may later update the symbol map contained in @p user_substitution_map
+     * by calling update_user_substitution_map().
+     *
+     * @param function A vector of symbolic expressions of type
+     * Differentiation::SD::Expression, representing the components of this
+     * Function.
+     *
+     * @param coordinate_symbols A tensor of symbols representing coordinates,
+     * used as input argument in the symbolic expressions contained in the
+     * @p function vector. The default @p coordinate_symbols is a
+     * Tensor<1,dim,Differentiation::SD::Expression>
+     * containing the symbols "x" for `dim` equal to one, "x", "y" for `dim`
+     * equal to two, and "x", "y", "z" for `dim` equal to three.
+     *
+     * @param time_symbol A symbolic variable representing time. It defaults
+     * to a symbolic variable named "t".
+     *
+     * @param user_substitution_map Any other symbol that may be contained in
+     * the symbolic function needs to be specified in this map. The map may be
+     * empty, and the functions may still contain unevaluated symbols, provided
+     * that you call update_user_substitution_map() and provide a replacement of
+     * all symbols except @p coordinate_symbols and @p time_symbol before any
+     * evaluation occurs.
+     */
+    SymbolicFunction(
+      const std::vector<Differentiation::SD::Expression> &function,
+      const Tensor<1, dim, Differentiation::SD::Expression>
+        &coordinate_symbols = get_default_coordinate_symbols(),
+      const Differentiation::SD::Expression &time_symbol =
+        Differentiation::SD::make_symbol("t"),
+      const Differentiation::SD::types::substitution_map
+        &user_substitution_map = {});
+
+    /**
+     * Constructor that takes a single string that describes the function
+     * expression as a semicolon separated list of expressions.
+     *
+     * The symbolic expression can use the default argument and the default
+     * symbolic time variable, plus any additional symbols that you may
+     * need, provided that you update the user substitution map that substitutes
+     * all of them before you try to evaluate the function or its derivatives,
+     * by calling update_user_substitution_map(), and that you provide all the
+     * additional function arguments of your function using the method
+     * set_additional_function_arguments().
+     */
+    SymbolicFunction(const std::string &expressions);
+
+    /**
+     * Store and apply the substitution map @p substitutions to each symbolic
+     * component of this Function object.
+     *
+     * Notice that this method will trigger a recomputation of the
+     * gradients, Hessians, and Laplacians of each component.
+     */
+    void
+    update_user_substitution_map(
+      const Differentiation::SD::types::substitution_map &substitutions);
+
+    /**
+     * Set the additional @p arguments to be substituted in next evaluation
+     * step.
+     *
+     * Notice that the @p arguments are substituted *after* evaluating the
+     * @p permanent_user_substitution_map, and after all derivatives are
+     * computed. If the additional arguments you pass still depend on the
+     * coordinate or time symbols, then evaluation of derivatives will result in
+     * a partial derivative evaluation.
+     *
+     * This method provides a way to evaluate functions that depend on more
+     * arguments than simply the coordinates and time. If you want to compute
+     * the total derivative w.r.t. to complicated symbolic expressions, you
+     * should call update_user_substitution_map() instead.
+     */
+    void
+    set_additional_function_arguments(
+      const Differentiation::SD::types::substitution_map &arguments);
+
+    /**
+     * Return a tensor of coordinate symbols that can be used to define the
+     * expressions of this symbolic function object.
+     *
+     * The default argument is a Tensor<1,dim,Differentiation::SD::Expression>
+     * containing the symbols "x" for `dim` equal to one, "x", "y" for `dim`
+     * equal to two, and "x", "y", "z" for `dim` equal to three.
+     */
+    static Tensor<1, dim, Differentiation::SD::Expression>
+    get_default_coordinate_symbols();
+
+    /**
+     * Get the actual arguments used for the coordinates in the symbolic
+     * function. This object does not include any user-defined arguments.
+     */
+    const Tensor<1, dim, Differentiation::SD::Expression> &
+    get_coordinate_symbols() const;
+
+    /**
+     * Get the actual symbolic time in use in this symbolic function.
+     */
+    const Differentiation::SD::Expression &
+    get_time_symbol() const;
+
+    /**
+     * Get the actual symbolic expressions used in this symbolic function.
+     */
+    const std::vector<Differentiation::SD::Expression> &
+    get_symbolic_function_expressions() const;
+
+    /**
+     * Get the currently stored @p user_substitution_map.
+     */
+    const Differentiation::SD::types::substitution_map &
+    get_user_substitution_map() const;
+
+    /**
+     * Return a SymbolicFunction object that represents the time derivative of
+     * this function. The spatial argument, the symbolic time, and the currently
+     * stored user substitution map are forwarded to the new function.
+     */
+    SymbolicFunction<dim, RangeNumberType>
+    time_derivative() const;
+
+    // documentation inherited from the base class
+    virtual RangeNumberType
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+
+    // documentation inherited from the base class
+    virtual Tensor<1, dim, RangeNumberType>
+    gradient(const Point<dim> & p,
+             const unsigned int component = 0) const override;
+
+    // documentation inherited from the base class
+    virtual RangeNumberType
+    laplacian(const Point<dim> & p,
+              const unsigned int component = 0) const override;
+
+    // documentation inherited from the base class
+    virtual SymmetricTensor<2, dim, RangeNumberType>
+    hessian(const Point<dim> & p,
+            const unsigned int component = 0) const override;
+
+    /**
+     * Print the stored arguments and function expression, as it would be
+     * evaluated when calling the method value().
+     */
+    template <typename StreamType>
+    StreamType &
+    print(StreamType &out) const;
+
+  private:
+    /**
+     * Return a substitution map that replaces the argument with the values of
+     * @p point, the symbolic time with the value of this->get_time(), and any
+     * additional arguments with the substitution map given by
+     * @p additional_function_arguments.
+     */
+    Differentiation::SD::types::substitution_map
+    create_evaluation_substitution_map(const Point<dim> &point) const;
+
+    /**
+     * Recompute the symbolic value of the function, applying the user
+     * substitution map. This may be an expensive computation, and it is called
+     * only if necessary.
+     */
+    void
+    update_values() const;
+
+    /**
+     * Recompute the symbolic gradient of the function, applying the user
+     * substitution map. This may be an expensive computation, and it is called
+     * only if necessary.
+     */
+    void
+    update_first_derivatives() const;
+
+    /**
+     * Recompute the symbolic Hessian and the symbolic Lapalacian of the
+     * function. This may be an expensive computation, and it is called
+     * only if necessary.
+     */
+    void
+    update_second_derivatives() const;
+
+    /**
+     * The components of this symbolic function, before any subustitution took
+     * place. This is immutable, and generated at construction time.
+     *
+     * Before any evaluation takes place, the @p user_substitution_map is
+     * applied to this object, and the result is stored in the internal variable
+     * function.
+     *
+     * During evaluation, the @p symbolic_coordinate, the @p symbolic_time, and
+     * any remaining symbols are substituted with the input evaluation point,
+     * the current time, and the content of @p additional_function_arguments.
+     */
+    const std::vector<Differentiation::SD::Expression> user_function;
+
+    /**
+     * Store the user substitution map used for expression substitutions. This
+     * may be updated with a call to update_user_substitution_map(). Notice that
+     * the function may still have unresolved symbols, provided that they are
+     * resolved by a call to set_additional_function_arguments().
+     */
+    Differentiation::SD::types::substitution_map user_substitution_map;
+
+    /**
+     * Store a user substitution map used for additional argument
+     * substitutions. This will be updated by a call to
+     * set_additional_function_arguments().
+     */
+    Differentiation::SD::types::substitution_map additional_function_arguments;
+
+    /**
+     * The actual components of this symbolic function. This is obtained from
+     * the @p user_function, after applying the @p user_substitution_map.
+     */
+    mutable std::vector<Differentiation::SD::Expression> function;
+
+    /**
+     * The gradients of each component of this symbolic function. This is
+     * obtained by computing the symbolic gradient of the object @p function,
+     * that is, after applying the @p user_substitution_map to @p user_function.
+     */
+    mutable std::vector<Tensor<1, dim, Differentiation::SD::Expression>>
+      function_gradient;
+
+    /**
+     * The Hessians of each component of this symbolic function. This is
+     * obtained by computing the symbolic Hessian of the object @p function,
+     * that is, after applying the @p user_substitution_map to @p user_function.
+     */
+    mutable std::vector<Tensor<2, dim, Differentiation::SD::Expression>>
+      function_hessian;
+
+    /**
+     * The Laplacians of each component of this symbolic function. This is
+     * obtained by computing the symbolic Laplacian of the object @p function,
+     * that is, after applying the @p user_substitution_map to @p user_function.
+     */
+    mutable std::vector<Differentiation::SD::Expression> function_laplacian;
+
+    /**
+     * The coordinate symbols argument of the function.
+     */
+    Tensor<1, dim, Differentiation::SD::Expression> coordinate_symbols;
+
+    /**
+     * The symbolic time argument of the function.
+     */
+    mutable Differentiation::SD::Expression time_symbol;
+  };
+
+  /**
+   * Allow output using the bitwise left shift operator.
+   */
+  template <int dim, typename RangeNumberType>
+  inline std::ostream &
+  operator<<(std::ostream &out, const SymbolicFunction<dim, RangeNumberType> &f)
+  {
+    return f.print(out);
+  }
+
+
+
+  // Inline and template functions
+  template <int dim, typename RangeNumberType>
+  template <typename StreamType>
+  StreamType &
+  SymbolicFunction<dim, RangeNumberType>::print(StreamType &out) const
+  {
+    for (unsigned int i = 0; i < dim; ++i)
+      out << coordinate_symbols[i] << ", ";
+    for (const auto it : additional_function_arguments)
+      out << it.first << ", ";
+    out << time_symbol << " -> " << user_function[0];
+    for (unsigned int i = 1; i < user_function.size(); ++i)
+      out << "; " << user_function[i];
+    if (!user_substitution_map.empty())
+      {
+        out << " # ( ";
+        std::string sep = "";
+        for (const auto it : user_substitution_map)
+          {
+            out << sep << it.first << " = " << it.second;
+            sep = ", ";
+          }
+        out << " )";
+      }
+    return out;
+  }
+#else
+  template <int dim, typename RangeNumberType = double>
+  class SymbolicFunction
+  {
+  public:
+    SymbolicFunction()
+    {
+      AssertThrow(
+        false,
+        ExcMessage(
+          "This class is not available if you did not enable SymEngine "
+          "when compiling deal.II."));
+    }
+  };
+#endif
+} // namespace Functions
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/base/symbolic_function.templates.h b/include/deal.II/base/symbolic_function.templates.h
new file mode 100644 (file)
index 0000000..6f30bd7
--- /dev/null
@@ -0,0 +1,301 @@
+// ---------------------------------------------------------------------
+//
+// 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_symbolic_function_templates_h
+#define dealii_symbolic_function_templates_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/symbolic_function.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Functions
+{
+#ifdef DEAL_II_WITH_SYMENGINE
+
+  namespace
+  {
+    inline std::vector<Differentiation::SD::Expression>
+    to_expressions(const std::vector<std::string> &strings)
+    {
+      std::vector<Differentiation::SD::Expression> ret(strings.size());
+      for (unsigned int i = 0; i < strings.size(); ++i)
+        ret[i] = Differentiation::SD::Expression(strings[i], true);
+      return ret;
+    }
+  } // namespace
+
+
+
+  template <int dim, typename RangeNumberType>
+  SymbolicFunction<dim, RangeNumberType>::SymbolicFunction(
+    const std::vector<Differentiation::SD::Expression> &   functions,
+    const Tensor<1, dim, Differentiation::SD::Expression> &argument,
+    const Differentiation::SD::Expression &                time,
+    const Differentiation::SD::types::substitution_map &user_substitution_map)
+    : Function<dim, RangeNumberType>(functions.size(), 0.0)
+    , user_function(functions)
+    , user_substitution_map(user_substitution_map)
+    , coordinate_symbols(argument)
+    , time_symbol(time)
+  {}
+
+
+
+  template <int dim, typename RangeNumberType>
+  SymbolicFunction<dim, RangeNumberType>::SymbolicFunction(
+    const std::string &expressions)
+    : Function<dim, RangeNumberType>(
+        Utilities::split_string_list(expressions, ';').size(),
+        0.0)
+    , user_function(
+        to_expressions(Utilities::split_string_list(expressions, ';')))
+    , coordinate_symbols(get_default_coordinate_symbols())
+    , time_symbol(Differentiation::SD::make_symbol("t"))
+  {}
+
+
+
+  template <int dim, typename RangeNumberType>
+  void
+  SymbolicFunction<dim, RangeNumberType>::update_values() const
+  {
+    // This is only necessary if the function vector is empty. A call to
+    // update_user_substitution_map() will clear all internal vectors.
+    if (function.empty())
+      {
+        function.resize(user_function.size());
+        for (unsigned int i = 0; i < user_function.size(); ++i)
+          {
+            function[i] = user_function[i].substitute(user_substitution_map);
+          }
+      }
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  void
+  SymbolicFunction<dim, RangeNumberType>::update_first_derivatives() const
+  {
+    // This is only necessary if the gradient vector is empty. A call to
+    // update_user_substitution_map() will clear all internal vectors.
+    if (function_gradient.empty())
+      {
+        update_values();
+        function_gradient.resize(user_function.size());
+        for (unsigned int i = 0; i < user_function.size(); ++i)
+          {
+            function_gradient[i] =
+              Differentiation::SD::differentiate(function[i],
+                                                 coordinate_symbols);
+          }
+      }
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  void
+  SymbolicFunction<dim, RangeNumberType>::update_second_derivatives() const
+  {
+    // This is only necessary if the gradient vector is empty. A call to
+    // update_user_substitution_map() will clear all internal vectors.
+    if (function_hessian.empty())
+      {
+        update_first_derivatives();
+        function_hessian.resize(user_function.size());
+        function_laplacian.resize(user_function.size());
+
+        for (unsigned int i = 0; i < user_function.size(); ++i)
+          {
+            function_hessian[i] =
+              Differentiation::SD::differentiate(function_gradient[i],
+                                                 coordinate_symbols);
+            function_laplacian[i] = trace(function_hessian[i]);
+          }
+      }
+  }
+
+
+  template <int dim, typename RangeNumberType>
+  void
+  SymbolicFunction<dim, RangeNumberType>::update_user_substitution_map(
+    const Differentiation::SD::types::substitution_map &substitutions)
+  {
+    user_substitution_map = substitutions;
+
+    // Now reset all internal vectors.
+    function.resize(0);
+    function_gradient.resize(0);
+    function_hessian.resize(0);
+    function_laplacian.resize(0);
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  void
+  SymbolicFunction<dim, RangeNumberType>::set_additional_function_arguments(
+    const Differentiation::SD::types::substitution_map &arguments)
+  {
+    additional_function_arguments = arguments;
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  Tensor<1, dim, Differentiation::SD::Expression>
+  SymbolicFunction<dim, RangeNumberType>::get_default_coordinate_symbols()
+  {
+    static_assert(dim <= 3, "Not implemented yet.");
+    const std::vector<std::string> names = {"x", "y", "z"};
+
+    Tensor<1, dim, Differentiation::SD::Expression> x;
+    for (unsigned int i = 0; i < dim; ++i)
+      x[i] = Differentiation::SD::make_symbol(names[i]);
+    return x;
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  const Tensor<1, dim, Differentiation::SD::Expression> &
+  SymbolicFunction<dim, RangeNumberType>::get_coordinate_symbols() const
+  {
+    return coordinate_symbols;
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  const Differentiation::SD::Expression &
+  SymbolicFunction<dim, RangeNumberType>::get_time_symbol() const
+  {
+    return time_symbol;
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  const std::vector<Differentiation::SD::Expression> &
+  SymbolicFunction<dim, RangeNumberType>::get_symbolic_function_expressions()
+    const
+  {
+    return user_function;
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  const Differentiation::SD::types::substitution_map &
+  SymbolicFunction<dim, RangeNumberType>::get_user_substitution_map() const
+  {
+    return user_substitution_map;
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  SymbolicFunction<dim, RangeNumberType>
+  SymbolicFunction<dim, RangeNumberType>::time_derivative() const
+  {
+    std::vector<Differentiation::SD::Expression> df_dt(user_function.size());
+    for (unsigned int i = 0; i < user_function.size(); ++i)
+      df_dt[i] = user_function[i].differentiate(time_symbol);
+    return SymbolicFunction<dim, RangeNumberType>(df_dt,
+                                                  coordinate_symbols,
+                                                  time_symbol,
+                                                  user_substitution_map);
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  RangeNumberType
+  SymbolicFunction<dim, RangeNumberType>::value(
+    const Point<dim> & p,
+    const unsigned int component) const
+  {
+    update_values();
+    AssertIndexRange(component, function.size());
+    return Differentiation::SD::substitute_and_evaluate<RangeNumberType>(
+      function[component], create_evaluation_substitution_map(p));
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  Tensor<1, dim, RangeNumberType>
+  SymbolicFunction<dim, RangeNumberType>::gradient(
+    const Point<dim> & p,
+    const unsigned int component) const
+  {
+    update_first_derivatives();
+    AssertIndexRange(component, function.size());
+    return Differentiation::SD::substitute_and_evaluate<RangeNumberType>(
+      function_gradient[component], create_evaluation_substitution_map(p));
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  RangeNumberType
+  SymbolicFunction<dim, RangeNumberType>::laplacian(
+    const Point<dim> & p,
+    const unsigned int component) const
+  {
+    update_second_derivatives();
+    AssertIndexRange(component, function.size());
+    return Differentiation::SD::substitute_and_evaluate<RangeNumberType>(
+      function_laplacian[component], create_evaluation_substitution_map(p));
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  SymmetricTensor<2, dim, RangeNumberType>
+  SymbolicFunction<dim, RangeNumberType>::hessian(
+    const Point<dim> & p,
+    const unsigned int component) const
+  {
+    update_second_derivatives();
+    AssertIndexRange(component, function.size());
+    return SymmetricTensor<2, dim, RangeNumberType>(
+      Differentiation::SD::substitute_and_evaluate<RangeNumberType>(
+        function_hessian[component], create_evaluation_substitution_map(p)));
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  Differentiation::SD::types::substitution_map
+  SymbolicFunction<dim, RangeNumberType>::create_evaluation_substitution_map(
+    const Point<dim> &p) const
+  {
+    auto map = additional_function_arguments;
+    Differentiation::SD::add_to_substitution_map(
+      map,
+      std::make_pair(time_symbol, this->get_time()),
+      std::make_pair(coordinate_symbols, p));
+    return map;
+  }
+#endif
+} // namespace Functions
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 2bc37ce8b4d8723476c5f2ad20d7fe69d7e1a0b1..be79de63f7b9ab51bb08c8d1cd81c7a443602dcd 100644 (file)
@@ -78,6 +78,7 @@ SET(_unity_include_src
   quadrature_selector.cc
   scalar_polynomials_base.cc
   subscriptor.cc
+  symbolic_function.cc
   table_handler.cc
   tensor_function.cc
   tensor_function_parser.cc
@@ -127,6 +128,7 @@ SET(_inst
   partitioner.cuda.inst.in
   polynomials_rannacher_turek.inst.in
   symmetric_tensor.inst.in
+  symbolic_function.inst.in
   tensor_function.inst.in
   tensor_function_parser.inst.in
   time_stepping.inst.in
diff --git a/source/base/symbolic_function.cc b/source/base/symbolic_function.cc
new file mode 100644 (file)
index 0000000..251050d
--- /dev/null
@@ -0,0 +1,26 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2015 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/symbolic_function.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+#ifdef DEAL_II_WITH_SYMENGINE
+namespace Functions
+{
+  // explicit instantiations
+#  include "symbolic_function.inst"
+} // namespace Functions
+#endif
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/base/symbolic_function.inst.in b/source/base/symbolic_function.inst.in
new file mode 100644 (file)
index 0000000..2afaa63
--- /dev/null
@@ -0,0 +1,25 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+for (S : REAL_SCALARS; dim : SPACE_DIMENSIONS)
+  {
+    template class SymbolicFunction<dim, S>;
+  }
+
+for (S : COMPLEX_SCALARS; dim : SPACE_DIMENSIONS)
+  {
+    template class SymbolicFunction<dim, S>;
+  }
diff --git a/tests/base/symbolic_function_01.cc b/tests/base/symbolic_function_01.cc
new file mode 100644 (file)
index 0000000..fd26ae6
--- /dev/null
@@ -0,0 +1,75 @@
+// ---------------------------------------------------------------------
+//
+// 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test the SymbolicFunction class
+
+#include <deal.II/base/symbolic_function.h>
+
+#include <map>
+
+#include "../tests.h"
+
+using namespace Differentiation::SD;
+
+template <int dim>
+Functions::SymbolicFunction<dim>
+test(const std::string &expression)
+{
+  auto x  = Functions::SymbolicFunction<dim>::get_default_coordinate_symbols();
+  auto t  = make_symbol("t");
+  auto f  = Expression(expression, true);
+  auto df = differentiate(f, x);
+  auto Hf = differentiate(df, x);
+
+  // Create a Function with a single component and an evaluation point
+  Functions::SymbolicFunction<dim> fun({f}, x);
+  Point<dim>                       p;
+  for (unsigned int i = 0; i < dim; ++i)
+    p[i] = i + 1.0;
+
+  // Output all symbolic stuff
+  deallog << "=========================================================="
+          << std::endl
+          << "dim = " << dim << ", Symengine" << std::endl
+          << "x: " << x << ", f: " << f << ", df: " << df << ", H:" << Hf
+          << std::endl;
+
+  // Output the function and its evaluation
+  deallog << "SymbolicFunction<dim>: " << fun << std::endl
+          << "p: " << p << ", f(p): " << fun.value(p)
+          << ", grad(f)(p): " << fun.gradient(p)
+          << ", laplacian(f)(p): " << fun.laplacian(p)
+          << ", H(f)(p): " << fun.hessian(p) << std::endl;
+
+  // Output its time derivative
+  auto fun_t = fun.time_derivative();
+  deallog << "Time derivative of SymbolicFunction<dim>: " << fun_t << std::endl
+          << "p: " << p << ", f_t(p): " << fun_t.value(p)
+          << ", grad(f_t)(p): " << fun_t.gradient(p)
+          << ", laplacian(f_t)(p): " << fun_t.laplacian(p)
+          << ", H(f_t)(p): " << fun_t.hessian(p) << std::endl;
+
+  return fun;
+}
+
+int
+main()
+{
+  initlog();
+
+  test<1>("x**2 + t*x");
+  test<2>("x+y+cos(exp(y*x*t))");
+  test<3>("tan(z*t**2) + atan(y/x)");
+}
diff --git a/tests/base/symbolic_function_01.with_symengine=true.output b/tests/base/symbolic_function_01.with_symengine=true.output
new file mode 100644 (file)
index 0000000..7523d0d
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL::==========================================================
+DEAL::dim = 1, Symengine
+DEAL::x: x, f: t*x + x**2, df: t + 2*x, H:2
+DEAL::SymbolicFunction<dim>: x, t -> t*x + x**2
+DEAL::p: 1.00000, f(p): 1.00000, grad(f)(p): 2.00000, laplacian(f)(p): 2.00000, H(f)(p): 2.00000
+DEAL::Time derivative of SymbolicFunction<dim>: x, t -> x
+DEAL::p: 1.00000, f_t(p): 1.00000, grad(f_t)(p): 1.00000, laplacian(f_t)(p): 0.00000, H(f_t)(p): 0.00000
+DEAL::==========================================================
+DEAL::dim = 2, Symengine
+DEAL::x: x y, f: x + y + cos(exp(t*x*y)), df: 1 - t*y*exp(t*x*y)*sin(exp(t*x*y)) 1 - t*x*exp(t*x*y)*sin(exp(t*x*y)), H:-t**2*y**2*exp(t*x*y)*sin(exp(t*x*y)) - t**2*y**2*exp(2*t*x*y)*cos(exp(t*x*y)) -t*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x*y*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x*y*exp(2*t*x*y)*cos(exp(t*x*y)) -t*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x*y*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x*y*exp(2*t*x*y)*cos(exp(t*x*y)) -t**2*x**2*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x**2*exp(2*t*x*y)*cos(exp(t*x*y))
+DEAL::SymbolicFunction<dim>: x, y, t -> x + y + cos(exp(t*x*y))
+DEAL::p: 1.00000 2.00000, f(p): 3.54030, grad(f)(p): 1.00000 1.00000, laplacian(f)(p): 0.00000, H(f)(p): 0.00000 0.00000 0.00000 0.00000
+DEAL::Time derivative of SymbolicFunction<dim>: x, y, t -> -x*y*exp(t*x*y)*sin(exp(t*x*y))
+DEAL::p: 1.00000 2.00000, f_t(p): -1.68294, grad(f_t)(p): -1.68294 -0.841471, laplacian(f_t)(p): 0.00000, H(f_t)(p): 0.00000 -0.841471 -0.841471 0.00000
+DEAL::==========================================================
+DEAL::dim = 3, Symengine
+DEAL::x: x y z, f: tan(t**2*z) + atan(y/x), df: -y/(x**2*(1 + y**2/x**2)) 1/(x*(1 + y**2/x**2)) t**2*(1 + tan(t**2*z)**2), H:-2*y**3/(x**5*(1 + y**2/x**2)**2) + 2*y/(x**3*(1 + y**2/x**2)) -1/(x**2*(1 + y**2/x**2)) + 2*y**2/(x**4*(1 + y**2/x**2)**2) 0 -1/(x**2*(1 + y**2/x**2)) + 2*y**2/(x**4*(1 + y**2/x**2)**2) -2*y/(x**3*(1 + y**2/x**2)**2) 0 0 0 2*t**4*tan(t**2*z)*(1 + tan(t**2*z)**2)
+DEAL::SymbolicFunction<dim>: x, y, z, t -> tan(t**2*z) + atan(y/x)
+DEAL::p: 1.00000 2.00000 3.00000, f(p): 1.10715, grad(f)(p): -0.400000 0.200000 0.00000, laplacian(f)(p): 2.77556e-17, H(f)(p): 0.160000 0.120000 0.00000 0.120000 -0.160000 0.00000 0.00000 0.00000 0.00000
+DEAL::Time derivative of SymbolicFunction<dim>: x, y, z, t -> 2*t*z*(1 + tan(t**2*z)**2)
+DEAL::p: 1.00000 2.00000 3.00000, f_t(p): 0.00000, grad(f_t)(p): 0.00000 0.00000 0.00000, laplacian(f_t)(p): 0.00000, H(f_t)(p): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
diff --git a/tests/base/symbolic_function_02.cc b/tests/base/symbolic_function_02.cc
new file mode 100644 (file)
index 0000000..3cfebae
--- /dev/null
@@ -0,0 +1,38 @@
+// ---------------------------------------------------------------------
+//
+// 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test the SymbolicFunction class: string constructor, and multiple components.
+
+#include <deal.II/base/symbolic_function.h>
+
+#include <map>
+
+#include "../tests.h"
+
+using namespace Differentiation::SD;
+
+int
+main()
+{
+  initlog();
+
+  Functions::SymbolicFunction<2> fun("x; y; t");
+  deallog << fun << std::endl;
+  Point<2> p(1, 2);
+  deallog << "p: " << p << std::endl
+          << "f_0(p): " << fun.value(p, 0) << std::endl
+          << "f_1(p): " << fun.value(p, 1) << std::endl
+          << "f_2(p): " << fun.value(p, 2) << std::endl;
+}
diff --git a/tests/base/symbolic_function_02.with_symengine=true.output b/tests/base/symbolic_function_02.with_symengine=true.output
new file mode 100644 (file)
index 0000000..7c07137
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::x, y, t -> x; y; t
+DEAL::p: 1.00000 2.00000
+DEAL::f_0(p): 1.00000
+DEAL::f_1(p): 2.00000
+DEAL::f_2(p): 0.00000
diff --git a/tests/base/symbolic_function_03.cc b/tests/base/symbolic_function_03.cc
new file mode 100644 (file)
index 0000000..76c2987
--- /dev/null
@@ -0,0 +1,50 @@
+// ---------------------------------------------------------------------
+//
+// 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test the SymbolicFunction class: string constructor, multiple components,
+// and updating constants
+
+#include <deal.II/base/symbolic_function.h>
+
+#include "../tests.h"
+
+using namespace Differentiation::SD;
+
+int
+main()
+{
+  initlog(0);
+
+  Functions::SymbolicFunction<2> fun("x; alpha*y; t*alpha^2");
+
+  auto sub = make_substitution_map("alpha", Expression("1.0", true));
+  fun.update_user_substitution_map(sub);
+
+  deallog << fun << std::endl;
+  Point<2> p(1, 2);
+  deallog << "p: " << p << std::endl
+          << "f_0(p): " << fun.value(p, 0) << std::endl
+          << "f_1(p): " << fun.value(p, 1) << std::endl
+          << "f_2(p): " << fun.value(p, 2) << std::endl;
+
+  sub = make_substitution_map(std::make_pair(make_symbol("alpha"), 2.0));
+  fun.update_user_substitution_map(sub);
+
+  deallog << fun << std::endl
+          << "p: " << p << std::endl
+          << "f_0(p): " << fun.value(p, 0) << std::endl
+          << "f_1(p): " << fun.value(p, 1) << std::endl
+          << "f_2(p): " << fun.value(p, 2) << std::endl;
+}
diff --git a/tests/base/symbolic_function_03.with_symengine=true.output b/tests/base/symbolic_function_03.with_symengine=true.output
new file mode 100644 (file)
index 0000000..d0f24a7
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::x, y, t -> x; y*alpha; t*alpha**2 # ( alpha = 1.0 )
+DEAL::p: 1.00000 2.00000
+DEAL::f_0(p): 1.00000
+DEAL::f_1(p): 2.00000
+DEAL::f_2(p): 0.00000
+DEAL::x, y, t -> x; y*alpha; t*alpha**2 # ( alpha = 2.0 )
+DEAL::p: 1.00000 2.00000
+DEAL::f_0(p): 1.00000
+DEAL::f_1(p): 4.00000
+DEAL::f_2(p): 0.00000
diff --git a/tests/base/symbolic_function_04.cc b/tests/base/symbolic_function_04.cc
new file mode 100644 (file)
index 0000000..0611d39
--- /dev/null
@@ -0,0 +1,53 @@
+// ---------------------------------------------------------------------
+//
+// 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test the SymbolicFunction class: string constructor, multiple components,
+// updating constants, and providing additional arguments.
+
+#include <deal.II/base/symbolic_function.h>
+
+#include "../tests.h"
+
+using namespace Differentiation::SD;
+
+int
+main()
+{
+  initlog(0);
+
+  Functions::SymbolicFunction<2> fun("x+beta; alpha*y; t*alpha^2");
+
+  auto sub = make_substitution_map("alpha", Expression("1.0", true));
+  fun.update_user_substitution_map(sub);
+
+  auto args = make_substitution_map(make_symbol("beta"), 0.0);
+  fun.set_additional_function_arguments(args);
+
+  deallog << fun << std::endl;
+  Point<2> p(1, 2);
+  deallog << "p: " << p << std::endl
+          << "f_0(p): " << fun.value(p, 0) << std::endl
+          << "f_1(p): " << fun.value(p, 1) << std::endl
+          << "f_2(p): " << fun.value(p, 2) << std::endl;
+
+  sub = make_substitution_map(std::make_pair(make_symbol("alpha"), 2.0));
+  fun.update_user_substitution_map(sub);
+
+  deallog << fun << std::endl
+          << "p: " << p << std::endl
+          << "f_0(p): " << fun.value(p, 0) << std::endl
+          << "f_1(p): " << fun.value(p, 1) << std::endl
+          << "f_2(p): " << fun.value(p, 2) << std::endl;
+}
diff --git a/tests/base/symbolic_function_04.with_symengine=true.output b/tests/base/symbolic_function_04.with_symengine=true.output
new file mode 100644 (file)
index 0000000..fcded2f
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::x, y, beta, t -> beta + x; y*alpha; t*alpha**2 # ( alpha = 1.0 )
+DEAL::p: 1.00000 2.00000
+DEAL::f_0(p): 1.00000
+DEAL::f_1(p): 2.00000
+DEAL::f_2(p): 0.00000
+DEAL::x, y, beta, t -> beta + x; y*alpha; t*alpha**2 # ( alpha = 2.0 )
+DEAL::p: 1.00000 2.00000
+DEAL::f_0(p): 1.00000
+DEAL::f_1(p): 4.00000
+DEAL::f_2(p): 0.00000

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.