--- /dev/null
+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)
* 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.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
quadrature_selector.cc
scalar_polynomials_base.cc
subscriptor.cc
+ symbolic_function.cc
table_handler.cc
tensor_function.cc
tensor_function_parser.cc
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>;
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.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)");
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.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;
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff