function_lib.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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_lib_h
#define dealii_function_lib_h


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

#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/table.h>

#include <array>

DEAL_II_NAMESPACE_OPEN

namespace Functions
{
  template <int dim>
  class SquareFunction : public Function<dim>
  {
  public:
    virtual double
    value(const Point<dim> &p, const unsigned int component = 0) const override;
    virtual void
    vector_value(const Point<dim> &p, Vector<double> &values) const override;
    virtual void
    value_list(const std::vector<Point<dim>> &points,
               std::vector<double> &          values,
               const unsigned int             component = 0) const override;
    virtual Tensor<1, dim>
    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>> &gradient) const override;
    virtual void
    gradient_list(const std::vector<Point<dim>> &points,
                  std::vector<Tensor<1, dim>> &  gradients,
                  const unsigned int             component = 0) const override;
    virtual double
    laplacian(const Point<dim> & p,
              const unsigned int component = 0) const override;
    virtual void
    laplacian_list(const std::vector<Point<dim>> &points,
                   std::vector<double> &          values,
                   const unsigned int             component = 0) const override;
  };



  template <int dim>
  class Q1WedgeFunction : public Function<dim>
  {
  public:
    virtual double
    value(const Point<dim> &p, const unsigned int component = 0) const override;

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

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

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

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

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

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

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



  template <int dim>
  class PillowFunction : public Function<dim>
  {
  public:
    PillowFunction(const double offset = 0.);

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

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

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

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

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

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

  private:
    const double offset;
  };



  template <int dim>
  class CosineFunction : public Function<dim>
  {
  public:
    CosineFunction(const unsigned int n_components = 1);

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

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

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

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

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

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

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

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

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



  template <int dim>
  class CosineGradFunction : public Function<dim>
  {
  public:
    CosineGradFunction();

    virtual double
    value(const Point<dim> &p, const unsigned int component) const override;
    virtual void
    vector_value(const Point<dim> &p, Vector<double> &values) const override;
    virtual void
    value_list(const std::vector<Point<dim>> &points,
               std::vector<double> &          values,
               const unsigned int             component) const override;

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

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

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

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

    virtual double
    laplacian(const Point<dim> &p, const unsigned int component) const override;
  };



  template <int dim>
  class ExpFunction : public Function<dim>
  {
  public:
    virtual double
    value(const Point<dim> &p, const unsigned int component = 0) const override;

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

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

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

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

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



  class LSingularityFunction : public Function<2>
  {
  public:
    virtual double
    value(const Point<2> &p, const unsigned int component = 0) const override;

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

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

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

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

    virtual void
    vector_gradient_list(
      const std::vector<Point<2>> &,
      std::vector<std::vector<Tensor<1, 2>>> &) const override;

    virtual double
    laplacian(const Point<2> &   p,
              const unsigned int component = 0) const override;

    virtual void
    laplacian_list(const std::vector<Point<2>> &points,
                   std::vector<double> &        values,
                   const unsigned int           component = 0) const override;
  };



  class LSingularityGradFunction : public Function<2>
  {
  public:
    LSingularityGradFunction();
    virtual double
    value(const Point<2> &p, const unsigned int component) const override;

    virtual void
    value_list(const std::vector<Point<2>> &points,
               std::vector<double> &        values,
               const unsigned int           component) const override;

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

    virtual Tensor<1, 2>
    gradient(const Point<2> &p, const unsigned int component) const override;

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

    virtual void
    vector_gradient_list(
      const std::vector<Point<2>> &,
      std::vector<std::vector<Tensor<1, 2>>> &) const override;

    virtual double
    laplacian(const Point<2> &p, const unsigned int component) const override;

    virtual void
    laplacian_list(const std::vector<Point<2>> &points,
                   std::vector<double> &        values,
                   const unsigned int           component) const override;
  };



  template <int dim>
  class SlitSingularityFunction : public Function<dim>
  {
  public:
    virtual double
    value(const Point<dim> &p, const unsigned int component = 0) const override;

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

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

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

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

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

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

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


  class SlitHyperSingularityFunction : public Function<2>
  {
  public:
    virtual double
    value(const Point<2> &p, const unsigned int component = 0) const override;

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

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

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

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

    virtual void
    vector_gradient_list(
      const std::vector<Point<2>> &,
      std::vector<std::vector<Tensor<1, 2>>> &) const override;

    virtual double
    laplacian(const Point<2> &   p,
              const unsigned int component = 0) const override;

    virtual void
    laplacian_list(const std::vector<Point<2>> &points,
                   std::vector<double> &        values,
                   const unsigned int           component = 0) const override;
  };



  template <int dim>
  class JumpFunction : public Function<dim>
  {
  public:
    JumpFunction(const Point<dim> &direction, const double steepness);

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

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

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

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

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

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

    std::size_t
    memory_consumption() const;

  protected:
    const Point<dim> direction;

    const double steepness;

    double angle;

    double sine;

    double cosine;
  };



  template <int dim>
  class FourierCosineFunction : public Function<dim>
  {
  public:
    FourierCosineFunction(const Tensor<1, dim> &fourier_coefficients);

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

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

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

  private:
    const Tensor<1, dim> fourier_coefficients;
  };



  template <int dim>
  class FourierSineFunction : public Function<dim>
  {
  public:
    FourierSineFunction(const Tensor<1, dim> &fourier_coefficients);

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

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

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

  private:
    const Tensor<1, dim> fourier_coefficients;
  };


  template <int dim>
  class FourierSineSum : public Function<dim>
  {
  public:
    FourierSineSum(const std::vector<Point<dim>> &fourier_coefficients,
                   const std::vector<double> &    weights);

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

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

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

  private:
    const std::vector<Point<dim>> fourier_coefficients;
    const std::vector<double>     weights;
  };



  template <int dim>
  class FourierCosineSum : public Function<dim>
  {
  public:
    FourierCosineSum(const std::vector<Point<dim>> &fourier_coefficients,
                     const std::vector<double> &    weights);

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

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

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

  private:
    const std::vector<Point<dim>> fourier_coefficients;
    const std::vector<double>     weights;
  };


  template <int dim>
  class CutOffFunctionBase : public Function<dim>
  {
  public:
    static const unsigned int no_component = numbers::invalid_unsigned_int;

    CutOffFunctionBase(
      const double       radius       = 1.,
      const Point<dim>   center       = Point<dim>(),
      const unsigned int n_components = 1,
      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
      const bool         integrate_to_one       = false,
      const double       unitary_integral_value = 1.0);

    virtual ~CutOffFunctionBase() = default;

    DEAL_II_DEPRECATED void
    new_center(const Point<dim> &p);

    DEAL_II_DEPRECATED void
    new_radius(const double r);

    virtual void
    set_center(const Point<dim> &p);

    virtual void
    set_radius(const double r);

    const Point<dim> &
    get_center() const;

    double
    get_radius() const;

    bool
    integrates_to_one() const;

  protected:
    Point<dim> center;

    double radius;

    const unsigned int selected;

    bool integrate_to_one;

    const double unitary_integral_value;

    double rescaling;
  };


  template <int dim>
  class CutOffFunctionTensorProduct : public CutOffFunctionBase<dim>
  {
  public:
    CutOffFunctionTensorProduct(
      double             radius       = 1.0,
      const Point<dim> & center       = Point<dim>(),
      const unsigned int n_components = 1,
      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
      const bool         integrate_to_one = false);

    template <template <int> class CutOffFunctionBaseType>
    void
    set_base();

    virtual void
    set_center(const Point<dim> &center) override;

    virtual void
    set_radius(const double radius) override;

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

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

  private:
    std::array<std::unique_ptr<CutOffFunctionBase<1>>, dim> base;

    bool initialized;
  };



  template <int dim>
  class CutOffFunctionLinfty : public CutOffFunctionBase<dim>
  {
  public:
    CutOffFunctionLinfty(
      const double radius             = 1.,
      const Point<dim>                = Point<dim>(),
      const unsigned int n_components = 1,
      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
      const bool         integrate_to_one = false);

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

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

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


  template <int dim>
  class CutOffFunctionW1 : public CutOffFunctionBase<dim>
  {
  public:
    CutOffFunctionW1(
      const double radius             = 1.,
      const Point<dim>                = Point<dim>(),
      const unsigned int n_components = 1,
      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
      const bool         integrate_to_one = false);

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

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

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


  template <int dim>
  class CutOffFunctionC1 : public CutOffFunctionBase<dim>
  {
  public:
    CutOffFunctionC1(
      const double radius             = 1.,
      const Point<dim>                = Point<dim>(),
      const unsigned int n_components = 1,
      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
      bool               integrate_to_one = false);

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

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

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

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


  template <int dim>
  class CutOffFunctionCinfty : public CutOffFunctionBase<dim>
  {
  public:
    CutOffFunctionCinfty(
      const double radius             = 1.,
      const Point<dim>                = Point<dim>(),
      const unsigned int n_components = 1,
      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
      bool               integrate_to_one = false);

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

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

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

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



  template <int dim>
  class Monomial : public Function<dim>
  {
  public:
    Monomial(const Tensor<1, dim> &exponents,
             const unsigned int    n_components = 1);

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

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

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

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

  private:
    const Tensor<1, dim> exponents;
  };



  template <int dim>
  class InterpolatedTensorProductGridData : public Function<dim>
  {
  public:
    InterpolatedTensorProductGridData(
      const std::array<std::vector<double>, dim> &coordinate_values,
      const Table<dim, double> &                  data_values);

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

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

  protected:
    TableIndices<dim>
    table_index_of_point(const Point<dim> &p) const;

    const std::array<std::vector<double>, dim> coordinate_values;

    const Table<dim, double> data_values;
  };


  template <int dim>
  class InterpolatedUniformGridData : public Function<dim>
  {
  public:
    InterpolatedUniformGridData(
      const std::array<std::pair<double, double>, dim> &interval_endpoints,
      const std::array<unsigned int, dim> &             n_subintervals,
      const Table<dim, double> &                        data_values);

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

  private:
    const std::array<std::pair<double, double>, dim> interval_endpoints;

    const std::array<unsigned int, dim> n_subintervals;

    const Table<dim, double> data_values;
  };


  template <int dim>
  class Polynomial : public Function<dim>
  {
  public:
    Polynomial(const Table<2, double> &   exponents,
               const std::vector<double> &coefficients);

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


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

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

  private:
    const Table<2, double> exponents;

    const std::vector<double> coefficients;
  };

#ifndef DOXYGEN



  // Template definitions
  template <int dim>
  template <template <int> class CutOffFunctionBaseType>
  void
  CutOffFunctionTensorProduct<dim>::set_base()
  {
    initialized = true;
    static_assert(
      std::is_base_of<CutOffFunctionBase<1>, CutOffFunctionBaseType<1>>::value,
      "You can only construct a CutOffFunctionTensorProduct function from "
      "a class derived from CutOffFunctionBase.");

    for (unsigned int i = 0; i < dim; ++i)
      base[i].reset(new CutOffFunctionBaseType<1>(this->radius,
                                                  Point<1>(this->center[i]),
                                                  this->n_components,
                                                  this->selected,
                                                  this->integrate_to_one));
  }



#endif

} // namespace Functions
DEAL_II_NAMESPACE_CLOSE

#endif