function.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 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_function_h
#define dealii_function_h


#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.h>
#include <deal.II/base/function_time.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 <functional>
#include <vector>

DEAL_II_NAMESPACE_OPEN


template <typename number>
class Vector;
template <int rank, int dim, typename Number>
class TensorFunction;

template <int dim, typename RangeNumberType = double>
class Function : public FunctionTime<
                   typename numbers::NumberTraits<RangeNumberType>::real_type>,
                 public Subscriptor
{
public:
  static const unsigned int dimension = dim;

  const unsigned int n_components;

  using time_type = typename FunctionTime<
    typename numbers::NumberTraits<RangeNumberType>::real_type>::time_type;

  Function(const unsigned int n_components = 1,
           const time_type    initial_time = 0.0);

  virtual ~Function() override = 0;

  Function &
  operator=(const Function &f);

  virtual RangeNumberType
  value(const Point<dim> &p, const unsigned int component = 0) const;

  virtual void
  vector_value(const Point<dim> &p, Vector<RangeNumberType> &values) const;

  virtual void
  value_list(const std::vector<Point<dim>> &points,
             std::vector<RangeNumberType> & values,
             const unsigned int             component = 0) const;

  virtual void
  vector_value_list(const std::vector<Point<dim>> &       points,
                    std::vector<Vector<RangeNumberType>> &values) const;

  virtual void
  vector_values(const std::vector<Point<dim>> &            points,
                std::vector<std::vector<RangeNumberType>> &values) const;

  virtual Tensor<1, dim, RangeNumberType>
  gradient(const Point<dim> &p, const unsigned int component = 0) const;

  virtual void
  vector_gradient(
    const Point<dim> &                            p,
    std::vector<Tensor<1, dim, RangeNumberType>> &gradients) const;

  virtual void
  gradient_list(const std::vector<Point<dim>> &               points,
                std::vector<Tensor<1, dim, RangeNumberType>> &gradients,
                const unsigned int component = 0) const;

  virtual void
  vector_gradients(
    const std::vector<Point<dim>> &                            points,
    std::vector<std::vector<Tensor<1, dim, RangeNumberType>>> &gradients) const;

  virtual void
  vector_gradient_list(
    const std::vector<Point<dim>> &                            points,
    std::vector<std::vector<Tensor<1, dim, RangeNumberType>>> &gradients) const;

  virtual RangeNumberType
  laplacian(const Point<dim> &p, const unsigned int component = 0) const;

  virtual void
  vector_laplacian(const Point<dim> &p, Vector<RangeNumberType> &values) const;

  virtual void
  laplacian_list(const std::vector<Point<dim>> &points,
                 std::vector<RangeNumberType> & values,
                 const unsigned int             component = 0) const;

  virtual void
  vector_laplacian_list(const std::vector<Point<dim>> &       points,
                        std::vector<Vector<RangeNumberType>> &values) const;

  virtual SymmetricTensor<2, dim, RangeNumberType>
  hessian(const Point<dim> &p, const unsigned int component = 0) const;

  virtual void
  vector_hessian(
    const Point<dim> &                                     p,
    std::vector<SymmetricTensor<2, dim, RangeNumberType>> &values) const;

  virtual void
  hessian_list(const std::vector<Point<dim>> &                        points,
               std::vector<SymmetricTensor<2, dim, RangeNumberType>> &values,
               const unsigned int component = 0) const;

  virtual void
  vector_hessian_list(
    const std::vector<Point<dim>> &                                     points,
    std::vector<std::vector<SymmetricTensor<2, dim, RangeNumberType>>> &values)
    const;


  std::size_t
  memory_consumption() const;
};


namespace Functions
{
  template <int dim, typename RangeNumberType = double>
  class ConstantFunction : public Function<dim, RangeNumberType>
  {
  public:
    ConstantFunction(const RangeNumberType value,
                     const unsigned int    n_components = 1);

    ConstantFunction(const std::vector<RangeNumberType> &values);

    ConstantFunction(const Vector<RangeNumberType> &values);

    ConstantFunction(const RangeNumberType *begin_ptr,
                     const unsigned int     n_components);

    virtual RangeNumberType
    value(const Point<dim> &p, const unsigned int component = 0) const override;

    virtual void
    vector_value(const Point<dim> &       p,
                 Vector<RangeNumberType> &return_value) const override;

    virtual void
    value_list(const std::vector<Point<dim>> &points,
               std::vector<RangeNumberType> & return_values,
               const unsigned int             component = 0) const override;

    virtual void
    vector_value_list(
      const std::vector<Point<dim>> &       points,
      std::vector<Vector<RangeNumberType>> &return_values) const override;

    virtual Tensor<1, dim, RangeNumberType>
    gradient(const Point<dim> & p,
             const unsigned int component = 0) const override;

    virtual void
    vector_gradient(
      const Point<dim> &                            p,
      std::vector<Tensor<1, dim, RangeNumberType>> &gradients) const override;

    virtual void
    gradient_list(const std::vector<Point<dim>> &               points,
                  std::vector<Tensor<1, dim, RangeNumberType>> &gradients,
                  const unsigned int component = 0) const override;

    virtual void
    vector_gradient_list(
      const std::vector<Point<dim>> &                            points,
      std::vector<std::vector<Tensor<1, dim, RangeNumberType>>> &gradients)
      const override;

    std::size_t
    memory_consumption() const;

  protected:
    std::vector<RangeNumberType> function_value_vector;
  };



  template <int dim, typename RangeNumberType = double>
  class ZeroFunction : public ConstantFunction<dim, RangeNumberType>
  {
  public:
    ZeroFunction(const unsigned int n_components = 1);
  };

} // namespace Functions

template <int dim, typename RangeNumberType = double>
using ConstantFunction DEAL_II_DEPRECATED =
  Functions::ConstantFunction<dim, RangeNumberType>;

template <int dim, typename RangeNumberType = double>
using ZeroFunction DEAL_II_DEPRECATED =
  Functions::ZeroFunction<dim, RangeNumberType>;



template <int dim, typename RangeNumberType = double>
class ComponentSelectFunction : public ConstantFunction<dim, RangeNumberType>
{
public:
  ComponentSelectFunction(const unsigned int    selected,
                          const RangeNumberType value,
                          const unsigned int    n_components);

  ComponentSelectFunction(const unsigned int selected,
                          const unsigned int n_components);

  ComponentSelectFunction(const std::pair<unsigned int, unsigned int> &selected,
                          const unsigned int n_components);


  virtual void
  substitute_function_value_with(
    const Functions::ConstantFunction<dim, RangeNumberType> &f);

  virtual void
  vector_value(const Point<dim> &       p,
               Vector<RangeNumberType> &return_value) const override;

  virtual void
  vector_value_list(
    const std::vector<Point<dim>> &       points,
    std::vector<Vector<RangeNumberType>> &values) const override;

  std::size_t
  memory_consumption() const;

protected:
  const std::pair<unsigned int, unsigned int> selected_components;
};



template <int dim, typename RangeNumberType = double>
class ScalarFunctionFromFunctionObject : public Function<dim, RangeNumberType>
{
public:
  ScalarFunctionFromFunctionObject(
    const std::function<RangeNumberType(const Point<dim> &)> &function_object);

  virtual RangeNumberType
  value(const Point<dim> &p, const unsigned int component = 0) const override;

private:
  const std::function<RangeNumberType(const Point<dim> &)> function_object;
};



template <int dim, typename RangeNumberType = double>
class VectorFunctionFromScalarFunctionObject
  : public Function<dim, RangeNumberType>
{
public:
  VectorFunctionFromScalarFunctionObject(
    const std::function<RangeNumberType(const Point<dim> &)> &function_object,
    const unsigned int selected_component,
    const unsigned int n_components);

  virtual RangeNumberType
  value(const Point<dim> &p, const unsigned int component = 0) const override;

  virtual void
  vector_value(const Point<dim> &       p,
               Vector<RangeNumberType> &values) const override;

private:
  const std::function<RangeNumberType(const Point<dim> &)> function_object;

  const unsigned int selected_component;
};


template <int dim, typename RangeNumberType = double>
class FunctionFromFunctionObjects : public Function<dim, RangeNumberType>
{
public:
  explicit FunctionFromFunctionObjects(const unsigned int n_components = 1,
                                       const double       initial_time = 0);

  FunctionFromFunctionObjects(
    const std::vector<std::function<RangeNumberType(const Point<dim> &)>>
      &          values,
    const double initial_time = 0.0);

  FunctionFromFunctionObjects(
    const std::vector<std::function<RangeNumberType(const Point<dim> &)>>
      &values,
    const std::vector<
      std::function<Tensor<1, dim, RangeNumberType>(const Point<dim> &)>>
      &          gradients,
    const double initial_time = 0.0);


  virtual RangeNumberType
  value(const Point<dim> &p, const unsigned int component = 0) const override;

  virtual Tensor<1, dim, RangeNumberType>
  gradient(const Point<dim> & p,
           const unsigned int component = 0) const override;

  void
  set_function_values(
    const std::vector<std::function<RangeNumberType(const Point<dim> &)>>
      &values);

  void
  set_function_gradients(
    const std::vector<
      std::function<Tensor<1, dim, RangeNumberType>(const Point<dim> &)>>
      &gradients);

private:
  std::vector<std::function<RangeNumberType(const Point<dim> &)>>
    function_values;

  std::vector<
    std::function<Tensor<1, dim, RangeNumberType>(const Point<dim> &)>>
    function_gradients;
};


template <int dim, typename RangeNumberType = double>
class VectorFunctionFromTensorFunction : public Function<dim, RangeNumberType>
{
public:
  VectorFunctionFromTensorFunction(
    const TensorFunction<1, dim, RangeNumberType> &tensor_function,
    const unsigned int                             selected_component = 0,
    const unsigned int                             n_components       = dim);

  virtual ~VectorFunctionFromTensorFunction() override = default;

  virtual RangeNumberType
  value(const Point<dim> &p, const unsigned int component = 0) const override;

  virtual void
  vector_value(const Point<dim> &       p,
               Vector<RangeNumberType> &values) const override;

  virtual void
  vector_value_list(
    const std::vector<Point<dim>> &       points,
    std::vector<Vector<RangeNumberType>> &value_list) const override;

private:
  const TensorFunction<1, dim, RangeNumberType> &tensor_function;

  const unsigned int selected_component;
};


#ifndef DOXYGEN
// icc 2018 complains about an undefined reference
// if we put this in the templates.h file
//
// The destructor is pure virtual so we can't default it
// in the declaration.
template <int dim, typename RangeNumberType>
inline Function<dim, RangeNumberType>::~Function() = default;
#endif


DEAL_II_NAMESPACE_CLOSE

#endif