symbolic_function.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2020 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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/enable_observer_pointer.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/function.h>
#include <deal.II/base/point.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
  template <int dim, typename RangeNumberType = double>
  class SymbolicFunction : public Function<dim, RangeNumberType>
  {
  public:
    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 = {});
 
    SymbolicFunction(const std::string &expressions);
 
    void
    update_user_substitution_map(
      const Differentiation::SD::types::substitution_map &substitutions);
 
    void
    set_additional_function_arguments(
      const Differentiation::SD::types::substitution_map &arguments);
 
    static Tensor<1, dim, Differentiation::SD::Expression>
    get_default_coordinate_symbols();
 
    const Tensor<1, dim, Differentiation::SD::Expression> &
    get_coordinate_symbols() const;
 
    const Differentiation::SD::Expression &
    get_time_symbol() const;
 
    const std::vector<Differentiation::SD::Expression> &
    get_symbolic_function_expressions() const;
 
    const Differentiation::SD::types::substitution_map &
    get_user_substitution_map() const;
 
    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;
 
    template <typename StreamType>
    StreamType &
    print(StreamType &out) const;
 
  private:
    Differentiation::SD::types::substitution_map
    create_evaluation_substitution_map(const Point<dim> &point) const;
 
    void
    update_values() const;
 
    void
    update_first_derivatives() const;
 
    void
    update_second_derivatives() const;
 
    const std::vector<Differentiation::SD::Expression> user_function;
 
    Differentiation::SD::types::substitution_map user_substitution_map;
 
    Differentiation::SD::types::substitution_map additional_function_arguments;
 
    mutable std::vector<Differentiation::SD::Expression> function;
 
    mutable std::vector<Tensor<1, dim, Differentiation::SD::Expression>>
      function_gradient;
 
    mutable std::vector<Tensor<2, dim, Differentiation::SD::Expression>>
      function_hessian;
 
    mutable std::vector<Differentiation::SD::Expression> function_laplacian;
 
    Tensor<1, dim, Differentiation::SD::Expression> coordinate_symbols;
 
    mutable Differentiation::SD::Expression time_symbol;
  };
 
  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 &argument_pair : additional_function_arguments)
      out << argument_pair.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 &substitution : user_substitution_map)
          {
            out << sep << substitution.first << " = " << substitution.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