common.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/analytic/static_wrapper.hpp>
#include <kernel/analytic/wrappers.hpp>
#include <kernel/util/math.hpp>
#include <kernel/geometry/cgal.hpp>
 
// includes, system
#include <initializer_list>
#include <deque>
#include <vector>
 
namespace FEAT
{
  namespace Analytic
  {
    namespace Common
    {
      template<typename Scalar_, typename DataType_>
      class TensorStatic
      {
      public:
        static DataType_ eval(DataType_ x)
        {
          return Scalar_::eval(x);
        }
 
        static DataType_ der_x(DataType_ x)
        {
          return Scalar_::der_x(x);
        }
 
        static DataType_ der_xx(DataType_ x)
        {
          return Scalar_::der_xx(x);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return Scalar_::eval(x) * Scalar_::eval(y);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y)
        {
          return Scalar_::der_x(x) * Scalar_::eval(y);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y)
        {
          return Scalar_::eval(x) * Scalar_::der_x(y);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::eval(x) * Scalar_::eval(y) * Scalar_::eval(z);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::der_x(x) * Scalar_::eval(y) * Scalar_::eval(z);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::eval(x) * Scalar_::der_x(y) * Scalar_::eval(z);
        }
 
        static DataType_ der_z(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::eval(x) * Scalar_::eval(y) * Scalar_::der_x(z);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y)
        {
          return Scalar_::der_xx(x) * Scalar_::eval(y);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y)
        {
          return Scalar_::eval(x) * Scalar_::der_xx(y);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y)
        {
          return Scalar_::der_x(x) * Scalar_::der_x(y);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y)
        {
          return Scalar_::der_x(x) * Scalar_::der_x(y);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::der_xx(x) * Scalar_::eval(y) * Scalar_::eval(z);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::eval(x) * Scalar_::der_xx(y) * Scalar_::eval(z);
        }
 
        static DataType_ der_zz(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::eval(x) * Scalar_::eval(y) * Scalar_::der_xx(z);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::der_x(x) * Scalar_::der_x(y) * Scalar_::eval(z);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::der_x(x) * Scalar_::der_x(y) * Scalar_::eval(z);
        }
 
        static DataType_ der_xz(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::der_x(x) * Scalar_::eval(y) * Scalar_::der_x(z);
        }
 
        static DataType_ der_zx(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::der_x(x) * Scalar_::eval(y) * Scalar_::der_x(z);
        }
 
        static DataType_ der_yz(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::eval(x) * Scalar_::der_x(y) * Scalar_::der_x(z);
        }
 
        static DataType_ der_zy(DataType_ x, DataType_ y, DataType_ z)
        {
          return Scalar_::eval(x) * Scalar_::der_x(y) * Scalar_::der_x(z);
        }
      }; // class TensorStatic<...>
 
      template<typename DataType_, int k_ = 1>
      class SineTensorStatic
      {
        static_assert(k_ > 0, "parameter k_ must be a positive integer");
 
      public:
        static DataType_ kpi()
        {
          return DataType_(k_) * Math::pi<DataType_>();
        }
 
        static DataType_ eval(DataType_ x)
        {
          return Math::sin(kpi() * x);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return Math::sin(kpi() * x) * Math::sin(kpi() * y);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y, DataType_ z)
        {
          return Math::sin(kpi() * x) * Math::sin(kpi() * y) * Math::sin(kpi() * z);
        }
 
        static DataType_ der_x(DataType_ x)
        {
          return kpi() * Math::cos(kpi() * x);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y)
        {
          return kpi() * Math::cos(kpi() * x) * Math::sin(kpi() * y);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y)
        {
          return kpi() * Math::sin(kpi() * x) * Math::cos(kpi() * y);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y, DataType_ z)
        {
          return kpi() * Math::cos(kpi() * x) * Math::sin(kpi() * y) * Math::sin(kpi() * z);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y, DataType_ z)
        {
          return kpi() * Math::sin(kpi() * x) * Math::cos(kpi() * y) * Math::sin(kpi() * z);
        }
 
        static DataType_ der_z(DataType_ x, DataType_ y, DataType_ z)
        {
          return kpi() * Math::sin(kpi() * x) * Math::sin(kpi() * y) * Math::cos(kpi() * z);
        }
 
        static DataType_ der_xx(DataType_ x)
        {
          return -Math::sqr(kpi()) * Math::sin(kpi() * x);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y)
        {
          return -Math::sqr(kpi()) * Math::sin(kpi() * x) * Math::sin(kpi() * y);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y)
        {
          return der_xx(x, y);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y)
        {
          return Math::sqr(kpi()) * Math::cos(kpi() * x) * Math::cos(kpi() * y);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y)
        {
          return der_xy(x, y);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y, DataType_ z)
        {
          return -Math::sqr(kpi()) * Math::sin(kpi() * x) * Math::sin(kpi() * y) * Math::sin(kpi() * z);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xx(x, y, z);
        }
 
        static DataType_ der_zz(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xx(x, y, z);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y, DataType_ z)
        {
          return Math::sqr(kpi()) * Math::cos(kpi() * x) * Math::cos(kpi() * y) * Math::sin(kpi() * z);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xy(x, y, z);
        }
 
        static DataType_ der_xz(DataType_ x, DataType_ y, DataType_ z)
        {
          return Math::sqr(kpi()) * Math::cos(kpi() * x) * Math::sin(kpi() * y) * Math::cos(kpi() * z);
        }
 
        static DataType_ der_zx(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xz(x, y, z);
        }
 
        static DataType_ der_yz(DataType_ x, DataType_ y, DataType_ z)
        {
          return Math::sqr(kpi()) * Math::sin(kpi() * x) * Math::cos(kpi() * y) * Math::cos(kpi() * z);
        }
 
        static DataType_ der_zy(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_yz(x, y, z);
        }
      }; // class SineTensorStatic<...>
 
      template<typename DataType_>
      using SineBubbleStatic = SineTensorStatic<DataType_, 1>;
 
      template<int dim_>
      using SineBubbleFunction = StaticWrapperFunction<dim_, SineBubbleStatic, true, true, true>;
 
      template<typename DataType_, int k_ = 1>
      class CosineTensorStatic
      {
        static_assert(k_ > 0, "parameter k_ must be a positive integer");
 
      public:
        static DataType_ kpi()
        {
          return DataType_(k_) * Math::pi<DataType_>();
        }
 
        static DataType_ eval(DataType_ x)
        {
          return Math::cos(kpi() * x);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return Math::cos(kpi() * x) * Math::cos(kpi() * y);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y, DataType_ z)
        {
          return Math::cos(kpi() * x) * Math::cos(kpi() * y) * Math::cos(kpi() * z);
        }
 
        static DataType_ der_x(DataType_ x)
        {
          return -kpi() * Math::sin(kpi() * x);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y)
        {
          return -kpi() * Math::sin(kpi() * x) * Math::cos(kpi() * y);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y)
        {
          return -kpi() * Math::cos(kpi() * x) * Math::sin(kpi() * y);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y, DataType_ z)
        {
          return -kpi() * Math::sin(kpi() * x) * Math::cos(kpi() * y) * Math::cos(kpi() * z);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y, DataType_ z)
        {
          return -kpi() * Math::cos(kpi() * x) * Math::sin(kpi() * y) * Math::cos(kpi() * z);
        }
 
        static DataType_ der_z(DataType_ x, DataType_ y, DataType_ z)
        {
          return -kpi() * Math::cos(kpi() * x) * Math::cos(kpi() * y) * Math::sin(kpi() * z);
        }
 
        static DataType_ der_xx(DataType_ x)
        {
          return -Math::sqr(kpi()) * Math::cos(kpi() * x);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y)
        {
          return -Math::sqr(kpi()) * Math::cos(kpi() * x) * Math::cos(kpi() * y);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y)
        {
          return der_xx(x, y);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y)
        {
          return Math::sqr(kpi()) * Math::sin(kpi() * x) * Math::sin(kpi() * y);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y)
        {
          return der_xy(x, y);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y, DataType_ z)
        {
          return -Math::sqr(kpi()) * Math::cos(kpi() * x) * Math::cos(kpi() * y) * Math::cos(kpi() * z);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xx(x, y, z);
        }
 
        static DataType_ der_zz(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xx(x, y, z);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y, DataType_ z)
        {
          return Math::sqr(kpi()) * Math::sin(kpi() * x) * Math::sin(kpi() * y) * Math::cos(kpi() * z);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xy(x, y, z);
        }
 
        static DataType_ der_xz(DataType_ x, DataType_ y, DataType_ z)
        {
          return Math::sqr(kpi()) * Math::sin(kpi() * x) * Math::cos(kpi() * y) * Math::sin(kpi() * z);
        }
 
        static DataType_ der_zx(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xz(x, y, z);
        }
 
        static DataType_ der_yz(DataType_ x, DataType_ y, DataType_ z)
        {
          return Math::sqr(kpi()) * Math::cos(kpi() * x) * Math::sin(kpi() * y) * Math::sin(kpi() * z);
        }
 
        static DataType_ der_zy(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_yz(x, y, z);
        }
      }; // class CosineTensorStatic<...>
 
      template<typename DataType_>
      using CosineWaveStatic = CosineTensorStatic<DataType_, 1>;
 
      template<int dim_>
      using CosineWaveFunction = StaticWrapperFunction<dim_, CosineWaveStatic, true, true, true>;
 
      template<typename DataType_>
      class ExpBubbleScalarStatic
      {
      public:
        static DataType_ eval(DataType_ x)
        {
          return (Math::exp(-Math::sqr(DataType_(2)*x - DataType_(1))) - Math::exp(-DataType_(1))) / (DataType_(1) - Math::exp(-DataType_(1)));
        }
 
        static DataType_ der_x(DataType_ x)
        {
          return (DataType_(8)*x - DataType_(4))*Math::exp(-Math::sqr(DataType_(2)*x - DataType_(1))) / (Math::exp(-DataType_(1)) - DataType_(1));
        }
 
        static DataType_ der_xx(DataType_ x)
        {
          return (DataType_(64)*x*(DataType_(1)-x)-DataType_(8)) * Math::exp(-Math::sqr(DataType_(2)*x - DataType_(1))) / (Math::exp(-DataType_(1)) - DataType_(1));
        }
      }; // class ExpBubbleScalarStatic<...>
 
      template<typename DataType_>
      using ExpBubbleStatic = TensorStatic<ExpBubbleScalarStatic<DataType_>, DataType_>;
 
      template<int dim_>
      using ExpBubbleFunction = StaticWrapperFunction<dim_, ExpBubbleStatic, true, true, true>;
 
      template<typename DataType_>
      class Q2BubbleScalarStatic
      {
      public:
        static DataType_ eval(DataType_ x)
        {
          return DataType_(4) * x * (DataType_(1) - x);
        }
 
        static DataType_ der_x(DataType_ x)
        {
          return DataType_(4) * (DataType_(1) - DataType_(2) * x);
        }
 
        static DataType_ der_xx(DataType_)
        {
          return -DataType_(8);
        }
      }; // class Q2BubbleScalarStatic<...>
 
      template<typename DataType_>
      using Q2BubbleStatic = TensorStatic<Q2BubbleScalarStatic<DataType_>, DataType_>;
 
      template<int dim_>
      using Q2BubbleFunction = StaticWrapperFunction<dim_, Q2BubbleStatic, true, true, true>;
 
      template<int dim_, typename DataType_ = Real>
      class ConstantFunction :
        public Analytic::Function
      {
      public:
        static constexpr int domain_dim = dim_;
        typedef Analytic::Image::Scalar ImageType;
 
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _value;
 
        public:
          explicit Evaluator(const ConstantFunction& function) :
            _value(DataType(function._value))
          {
          }
 
          ValueType value(const PointType& DOXY(point))
          {
            return _value;
          }
 
          GradientType gradient(const PointType& DOXY(point))
          {
            return GradientType::null();
          }
 
          HessianType hessian(const PointType& DOXY(point))
          {
            return HessianType::null();
          }
        }; // class ConstantFunction::Evaluator<...>
 
      private:
        DataType_ _value;
 
      public:
        explicit ConstantFunction(DataType_ value = DataType_(0)) :
          _value(value)
        {
        }
      }; // class ConstantFunction
 
      template<int dim_, typename DataType_ = Real>
      class ConstantVectorFunction :
        public Analytic::Function
      {
      public:
        static constexpr int domain_dim = dim_;
        typedef Analytic::Image::Vector<domain_dim> ImageType;
 
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const ValueType _value_vec;
 
        public:
          explicit Evaluator(const ConstantVectorFunction& function) :
            _value_vec(function._value_vector)
          {
          }
 
          ValueType value(const PointType& DOXY(point))
          {
            return _value_vec;
          }
 
          GradientType gradient(const PointType& DOXY(point))
          {
            return GradientType::null();
          }
 
          HessianType hessian(const PointType& DOXY(point))
          {
            return HessianType::null();
          }
        }; // class ConstantVectorFunction::Evaluator<...>
 
      private:
        Tiny::Vector<DataType_, domain_dim> _value_vector;
 
      public:
        explicit ConstantVectorFunction(DataType_ value = DataType_(0)) :
          _value_vector(Tiny::Vector<DataType_, domain_dim>(value))
        {
        }
 
        explicit ConstantVectorFunction(Tiny::Vector<DataType_, domain_dim> vec) :
          _value_vector(vec)
        {
        }
 
        explicit ConstantVectorFunction(DataType_ value_x, DataType_ value_y) :
          _value_vector()
        {
          _value_vector[0] = value_x;
          if constexpr(dim_ > 1)
            _value_vector[1] = value_y;
        }
 
        explicit ConstantVectorFunction(DataType_ value_x, DataType_ value_y, DataType_ value_z) :
          _value_vector()
        {
          _value_vector[0] = value_x;
          if constexpr(dim_ > 1)
            _value_vector[1] = value_y;
          if constexpr(dim_ > 2)
            _value_vector[2] = value_z;
        }
      }; // class ConstantVectorFunction
 
      template<typename AnalyticFunctionType1, typename AnalyticFunctionType2>
      class MinOfTwoFunctions :
        public Analytic::Function
      {
      public:
        static_assert(AnalyticFunctionType1::domain_dim == AnalyticFunctionType2::domain_dim, "domain dimension mismatch");
 
        static_assert(std::is_same<typename AnalyticFunctionType1::ImageType, Analytic::Image::Scalar>::value, "invalid image type");
        static_assert(std::is_same<typename AnalyticFunctionType2::ImageType, Analytic::Image::Scalar>::value, "invalid image type");
 
        static constexpr int domain_dim = AnalyticFunctionType1::domain_dim;
        typedef Analytic::Image::Scalar ImageType;
 
        static constexpr bool can_value = (AnalyticFunctionType1::can_value && AnalyticFunctionType2::can_value);
        static constexpr bool can_grad = (AnalyticFunctionType1::can_grad && AnalyticFunctionType2::can_grad);
        static constexpr bool can_hess = (AnalyticFunctionType1::can_hess && AnalyticFunctionType2::can_hess);
 
      public:
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          typename AnalyticFunctionType1::template Evaluator<EvalTraits_> _f1_eval;
          typename AnalyticFunctionType2::template Evaluator<EvalTraits_> _f2_eval;
 
          const DataType _eps;
 
        public:
          explicit Evaluator(const MinOfTwoFunctions& function) :
            _f1_eval(function._f1),
            _f2_eval(function._f2),
            _eps(Math::eps<DataType>())
          {
          }
 
          ValueType value(const PointType& point)
          {
            return Math::min(_f1_eval.value(point), _f2_eval.value(point));
          }
 
          GradientType gradient(const PointType& point)
          {
            ValueType val1 = _f1_eval.value(point);
            ValueType val2 = _f2_eval.value(point);
 
            if(Math::abs(val1-val2) < _eps)
              return GradientType::null();
            else if(val1 < val2)
              return _f1_eval.gradient(point);
            else
              return _f2_eval.gradient(point);
          }
 
          HessianType hessian(const PointType& point)
          {
            ValueType val1 = _f1_eval.value(point);
            ValueType val2 = _f2_eval.value(point);
 
            if(Math::abs(val1-val2) < _eps)
              return HessianType::null();
            else if(val1 < val2)
              return _f1_eval.hessian(point);
            else
              return _f2_eval.hessian(point);
          }
        }; // class MinOfTwoFunctions::Evaluator<...>
 
      private:
        const AnalyticFunctionType1& _f1;
        const AnalyticFunctionType2& _f2;
 
      public:
        explicit MinOfTwoFunctions(const AnalyticFunctionType1& f1_, const AnalyticFunctionType2& f2_) :
          _f1(f1_),
          _f2(f2_)
        {
        }
      }; // class MinOfTwoFunctions
 
 
      template<typename DataType_>
      class HeavisideStatic
      {
      public:
 
        static DataType_ eval(DataType_ x)
        {
          return ( x < DataType_(0) ) ? DataType_(0) : DataType_(1);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return eval(x)*eval(y);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y, DataType_ z)
        {
          return eval(x)*eval(y)*eval(z);
        }
 
      }; // class HeavisideStatic<...>
 
      template<int dim_>
      using HeavisideFunction = StaticWrapperFunction<dim_, HeavisideStatic, true, false, false>;
 
 
      template<typename DataType_>
      class HeavisideRegStatic
      {
      public:
 
        static DataType_ eval(DataType_ x)
        {
          return ( x < DataType_(0) ) ? DataType_(0) : DataType_(2)*(Math::cosh(x) - DataType_(1));
        }
 
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return eval(x)*eval(y);
        }
 
        static DataType_ eval(DataType_ x, DataType_ y, DataType_ z)
        {
          return eval(x)*eval(y)*eval(z);
        }
 
        static DataType_ der_x(DataType_ x)
        {
          return ( x < DataType_(0) ) ? DataType_(0) : DataType_(2)*Math::sinh(x);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y)
        {
          return der_x(x)*eval(y);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y)
        {
          return eval(x)*der_x(y);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_x(x)*eval(y)*eval(z);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y, DataType_ z)
        {
          return eval(x)*der_x(y)*eval(z);
        }
 
        static DataType_ der_z(DataType_ x, DataType_ y, DataType_ z)
        {
          return eval(x)*eval(y)*der_x(z);
        }
 
        static DataType_ der_xx(DataType_ x)
        {
          return ( x < DataType_(0) ) ? DataType_(0) : DataType_(2)*Math::cosh(x);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y)
        {
          return der_xx(x)*eval(y);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y)
        {
          return eval(x)*der_xx(y);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y)
        {
          return der_x(x)*der_x(y);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y)
        {
          return der_xy(x, y);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xx(x)*eval(y)*eval(z);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y, DataType_ z)
        {
          return eval(x)*der_xx(y)*eval(z);
        }
 
        static DataType_ der_zz(DataType_ x, DataType_ y, DataType_ z)
        {
          return eval(x)*eval(y)*der_xx(z);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_x(x)*der_x(y)*eval(z);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_x(x)*der_x(y)*eval(z);
        }
 
        static DataType_ der_xz(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_x(x)*eval(y)*der_x(z);
        }
 
        static DataType_ der_zx(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_xz(x, y, z);
        }
 
        static DataType_ der_yz(DataType_ x, DataType_ y, DataType_ z)
        {
          return eval(x)*der_x(y)*der_x(z);
        }
 
        static DataType_ der_zy(DataType_ x, DataType_ y, DataType_ z)
        {
          return der_yz(x, y, z);
        }
 
      }; // class HeavisideRegStatic<...>
 
      template<int dim_>
      using HeavisideRegFunction = StaticWrapperFunction<dim_, HeavisideRegStatic, true, false, false>;
 
      template<typename DataType_>
      class PolynomialFunction1D :
        public Analytic::Function
      {
      public:
        static constexpr int domain_dim = 1;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          std::vector<DataType> _coeff;
 
        public:
          explicit Evaluator(const PolynomialFunction1D& function)
          {
            for(auto it = function._coeff.begin(); it != function._coeff.end(); ++it)
              _coeff.push_back(DataType(*it));
          }
 
          ValueType value(const PointType& point) const
          {
            // evaluate polynomial via horner scheme
            DataType x = point[0];
            DataType y = DataType(0);
            for(std::size_t k(_coeff.size()); k > std::size_t(0); )
              y = x * y + _coeff[--k];
 
            return y;
          }
 
          GradientType gradient(const PointType& point) const
          {
            std::size_t k = _coeff.size();
            if(k <= std::size_t(0))
              return GradientType::null();
 
            // evaluate polynomial via horner scheme
            DataType x = point[0];
            DataType y = DataType(0);
            for( ; (--k) > std::size_t(0); )
              y = x * y + (_coeff[k] * DataType(k));
 
            return GradientType(y);
          }
 
          HessianType hessian(const PointType& point) const
          {
            std::size_t k = _coeff.size();
            if(k <= std::size_t(1))
              return HessianType::null();
 
            // evaluate polynomial via horner scheme
            DataType x = point[0];
            DataType y = DataType(0);
            for( ; (--k) > std::size_t(1); )
              y = x * y + (_coeff[k] * DataType(k*(k-1)));
 
            return HessianType(y);
          }
        }; // class PolynomialFunction1D::Evaluator
 
      private:
        std::vector<DataType_> _coeff;
 
      public:
        PolynomialFunction1D() :
          _coeff()
        {
        }
 
        explicit PolynomialFunction1D(std::initializer_list<DataType_> coeff) :
          _coeff(coeff)
        {
        }
 
        Index set_coeff(Index degree, DataType_ coeff)
        {
          // check degree
          if(Index(_coeff.size()) <= degree)
            _coeff.resize(std::size_t(degree+1), DataType_(0));
 
          // set coefficient
          _coeff.at(std::size_t(degree)) = coeff;
 
          // return degree+1
          return Index(_coeff.size());
        }
      }; // class PolynomialFunction1D
 
 
      template<typename DataType_>
      class BazaraaShettyStatic
      {
      public:
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return Math::sqr(Math::sqr(x - DataType_(2))) + Math::sqr(x - DataType_(2)*y);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y)
        {
          return DataType_(4)*(x - DataType_(2))*Math::sqr(x - DataType_(2)) + DataType_(2)*(x - DataType_(2)*y);
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y)
        {
          return DataType_(4)*(DataType_(2)*y - x);
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ DOXY(y))
        {
          return DataType_(12)*Math::sqr(x - DataType_(2)) + DataType_(2);
        }
 
        static DataType_ der_yy(DataType_ DOXY(x), DataType_ DOXY(y))
        {
          return DataType_(8);
        }
 
        static DataType_ der_xy(DataType_ DOXY(x), DataType_ DOXY(y))
        {
          return DataType_(-4);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y)
        {
          return der_xy(x,y);
        }
 
      }; // class BazaraaShettyStatic<...>
 
 
      template<typename DataType_>
      class GoldsteinPriceStatic
      {
      public:
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return (DataType_(1) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19))) * (DataType_(30) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18)));
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y)
        {
          return (DataType_(2) * (DataType_(1) + x + y) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19)) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(6) * x + DataType_(6) * y - DataType_(14))) * (DataType_(30) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18))) + (DataType_(1) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19))) * (DataType_(4) * (DataType_(2) * x - DataType_(3) * y) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18)) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(24) * x - DataType_(36) * y - DataType_(32)));
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y)
        {
          return (DataType_(2) * (DataType_(1) + x + y) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19)) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(6) * x + DataType_(6) * y - DataType_(14))) * (DataType_(30) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18))) + (DataType_(1) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19))) * (-DataType_(6) * (DataType_(2) * x - DataType_(3) * y) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18)) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (-DataType_(36) * x + DataType_(54) * y + DataType_(48)));
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y)
        {
          return DataType_((DataType_(6) * Math::pow(x, DataType_(2)) + DataType_(12) * x * y + DataType_(6) * Math::pow(y, DataType_(2)) - DataType_(28) * x - DataType_(28) * y + DataType_(38) + DataType_(4) * (DataType_(1) + x + y) * (DataType_(6) * x + DataType_(6) * y - DataType_(14)) + DataType_(6) * Math::pow(DataType_(1) + x + y, DataType_(2))) * (DataType_(30) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18))) + DataType_(2) * (DataType_(2) * (DataType_(1) + x + y) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19)) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(6) * x + DataType_(6) * y - DataType_(14))) * (DataType_(4) * (DataType_(2) * x - DataType_(3) * y) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18)) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(24) * x - DataType_(36) * y - DataType_(32))) + (DataType_(1) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19))) * (DataType_(96) * Math::pow(x, DataType_(2)) - DataType_(288) * x * y + DataType_(216) * Math::pow(y, DataType_(2)) - DataType_(256) * x + DataType_(384) * y + DataType_(144) + DataType_(8) * (DataType_(2) * x - DataType_(3) * y) * (DataType_(24) * x - DataType_(36) * y - DataType_(32)) + DataType_(24) * Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2))));
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y)
        {
          return (DataType_(6) * Math::pow(x, DataType_(2)) + DataType_(12) * x * y + DataType_(6) * Math::pow(y, DataType_(2)) - DataType_(28) * x - DataType_(28) * y + DataType_(38) + DataType_(4) * (DataType_(1) + x + y) * (DataType_(6) * x + DataType_(6) * y - DataType_(14)) + DataType_(6) * Math::pow(DataType_(1) + x + y, DataType_(2))) * (DataType_(30) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18))) + DataType_(2) * (DataType_(2) * (DataType_(1) + x + y) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19)) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(6) * x + DataType_(6) * y - DataType_(14))) * (-DataType_(6) * (DataType_(2) * x - DataType_(3) * y) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18)) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (-DataType_(36) * x + DataType_(54) * y + DataType_(48))) + (DataType_(1) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19))) * (DataType_(216) * Math::pow(x, DataType_(2)) - DataType_(648) * x * y + DataType_(486) * Math::pow(y, DataType_(2)) - DataType_(576) * x + DataType_(864) * y + DataType_(324) - DataType_(12) * (DataType_(2) * x - DataType_(3) * y) * (-DataType_(36) * x + DataType_(54) * y + DataType_(48)) + DataType_(54) * Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)));
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y)
        {
          return (DataType_(6) * Math::pow(x, DataType_(2)) + DataType_(12) * x * y + DataType_(6) * Math::pow(y, DataType_(2)) - DataType_(28) * x - DataType_(28) * y + DataType_(38) + DataType_(4) * (DataType_(1) + x + y) * (DataType_(6) * x + DataType_(6) * y - DataType_(14)) + DataType_(6) * Math::pow(DataType_(1) + x + y, DataType_(2))) * (DataType_(30) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18))) + (DataType_(2) * (DataType_(1) + x + y) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19)) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(6) * x + DataType_(6) * y - DataType_(14))) * (DataType_(4) * (DataType_(2) * x - DataType_(3) * y) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18)) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (DataType_(24) * x - DataType_(36) * y - DataType_(32))) + (DataType_(2) * (DataType_(1) + x + y) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19)) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(6) * x + DataType_(6) * y - DataType_(14))) * (-DataType_(6) * (DataType_(2) * x - DataType_(3) * y) * (DataType_(12) * Math::pow(x, DataType_(2)) - DataType_(36) * x * y + DataType_(27) * Math::pow(y, DataType_(2)) - DataType_(32) * x + DataType_(48) * y + DataType_(18)) + Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)) * (-DataType_(36) * x + DataType_(54) * y + DataType_(48))) + (DataType_(1) + Math::pow(DataType_(1) + x + y, DataType_(2)) * (DataType_(3) * Math::pow(x, DataType_(2)) + DataType_(6) * x * y + DataType_(3) * Math::pow(y, DataType_(2)) - DataType_(14) * x - DataType_(14) * y + DataType_(19))) * (-DataType_(144) * Math::pow(x, DataType_(2)) + DataType_(432) * x * y - DataType_(324) * Math::pow(y, DataType_(2)) + DataType_(384) * x - DataType_(576) * y - DataType_(216) - DataType_(6) * (DataType_(2) * x - DataType_(3) * y) * (DataType_(24) * x - DataType_(36) * y - DataType_(32)) + DataType_(4) * (DataType_(2) * x - DataType_(3) * y) * (-DataType_(36) * x + DataType_(54) * y + DataType_(48)) - DataType_(36) * Math::pow(DataType_(2) * x - DataType_(3) * y, DataType_(2)));
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y)
        {
          return der_xy(x, y);
        }
 
      }; // class GoldsteinPriceStatic<...>
 
      using GoldsteinPriceFunction = StaticWrapperFunction<2, GoldsteinPriceStatic, true, true, true>;
 
      using BazaraaShettyFunction = StaticWrapperFunction<2, BazaraaShettyStatic, true, true, true>;
 
 
      template<typename DataType_>
      class HimmelblauStatic
      {
      public:
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return Math::sqr(Math::sqr(x) + y - DataType_(11))
            + Math::sqr( x + Math::sqr(y) - DataType_(7));
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y)
        {
          return DataType_(4)*x*(Math::sqr(x) + y - DataType_(11)) + DataType_(2)*(x + Math::sqr(y) - DataType_(7));
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y)
        {
          return DataType_(2)* (Math::sqr(x) + y - DataType_(11)) + DataType_(4)*y*(x + Math::sqr(y) - DataType_(7));
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y)
        {
          return DataType_(12)*Math::sqr(x) + DataType_(4)*y - DataType_(42);
        }
 
        static DataType_ der_yy(DataType_ x, DataType_ y)
        {
          return DataType_(4)*x + DataType_(12)*Math::sqr(y) - DataType_(26);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ y)
        {
          return DataType_(4)*(x + y);
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y)
        {
          return der_xy(x,y);
        }
 
      }; // class HimmelblauStatic<...>
 
      using HimmelblauFunction = StaticWrapperFunction<2, HimmelblauStatic, true, true, true>;
 
      template<typename DataType_>
      class RosenbrockStatic
      {
      public:
        static DataType_ eval(DataType_ x, DataType_ y)
        {
          return DataType_(100)*Math::sqr( y - Math::sqr(x)) + Math::sqr(DataType_(1) - x);
        }
 
        static DataType_ der_x(DataType_ x, DataType_ y)
        {
          return -DataType_(400)*x*(y - Math::sqr(x) ) + DataType_(2)*(x - DataType_(1));
        }
 
        static DataType_ der_y(DataType_ x, DataType_ y)
        {
          return DataType_(200)*(y - Math::sqr(x));
        }
 
        static DataType_ der_xx(DataType_ x, DataType_ y)
        {
          return DataType_(1200)*Math::sqr(x) - DataType_(400)*y + DataType_(2);
        }
 
        static DataType_ der_yy(DataType_ DOXY(x), DataType_ DOXY(y))
        {
          return DataType_(200);
        }
 
        static DataType_ der_xy(DataType_ x, DataType_ DOXY(y))
        {
          return -DataType_(400)*x;
        }
 
        static DataType_ der_yx(DataType_ x, DataType_ y)
        {
          return der_xy(x,y);
        }
 
      }; // class RosenbrockStatic<...>
 
      using RosenbrockFunction = StaticWrapperFunction<2, RosenbrockStatic, true, true, true>;
 
      template<typename DataType_>
      class ParProfileBase :
        public Analytic::Function
      {
      public:
        static constexpr int domain_dim = 2;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr bool can_value = true;
 
      protected:
        // coordinates of line segment
        DataType_ _x0, _y0, _x1, _y1;
        // maximum value
        DataType_ _vmax;
 
      public:
        ParProfileBase() :
          _x0(0.0), _y0(0.0), _x1(0.0), _y1(1.0), _vmax(1.0)
        {
        }
 
        explicit ParProfileBase(DataType_ x0, DataType_ y0, DataType_ x1, DataType_ y1, DataType_ vmax = DataType_(1.0)) :
          _x0(x0), _y0(y0), _x1(x1), _y1(y1), _vmax(vmax)
        {
        }
 
        bool parse(const String& sbc)
        {
          auto li = sbc.find_first_of('(');
          auto ri = sbc.find_last_of(')');
          if((li == sbc.npos) || (ri == sbc.npos))
            return false;
 
          std::deque<String> sv = sbc.substr(li+1, ri-li-1).split_by_string(",");
          if((sv.size() < std::size_t(2)) || (sv.size() > std::size_t(3)))
            return false;
 
          std::deque<String> sv0 = sv[0].trim().split_by_whitespaces();
          std::deque<String> sv1 = sv[1].trim().split_by_whitespaces();
          if(sv0.size() != std::size_t(2)) return false;
          if(sv1.size() != std::size_t(2)) return false;
 
          if(!sv0[0].parse(_x0)) return false;
          if(!sv0[1].parse(_y0)) return false;
          if(!sv1[0].parse(_x1)) return false;
          if(!sv1[1].parse(_y1)) return false;
 
          if((sv.size() == std::size_t(3)) && !sv.back().parse(_vmax))
            return false;
 
          return true;
        }
      }; // class ParProfileBase
 
      template<typename DataType_>
      class ParProfileScalar :
        public ParProfileBase<DataType_>
      {
      public:
        static constexpr int domain_dim = 2;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
        using ParProfileBase<DataType_>::ParProfileBase;
 
        template<typename Traits_>
        class Evaluator :
          public Analytic::Function::Evaluator<Traits_>
        {
        protected:
          typedef typename Traits_::DataType DataType;
          typedef typename Traits_::PointType PointType;
          typedef typename Traits_::ValueType ValueType;
          typedef typename Traits_::GradientType GradientType;
          typedef typename Traits_::HessianType HessianType;
 
          PointType _vo, _ve;
          DataType _den, _vmax;
 
        public:
          explicit Evaluator(const ParProfileScalar& function)
          {
            _vo[0] = DataType(function._x0);
            _vo[1] = DataType(function._y0);
            _ve[0] = DataType(function._x1 - function._x0);
            _ve[1] = DataType(function._y1 - function._y0);
            _den = DataType(1) / Tiny::dot(_ve, _ve);
            _vmax = DataType(function._vmax);
          }
 
          ValueType value(const PointType& point)
          {
            // project point onto line segment
            const DataType x = Math::clamp(Tiny::dot(point - _vo, _ve) * _den, DataType(0), DataType(1));
            // compute function value
            return _vmax * DataType(4) * x * (DataType(1) - x);
          }
 
          GradientType gradient(const PointType& point)
          {
            // project point onto line segment
            const DataType x = Tiny::dot(point - _vo, _ve) * _den;
 
            // Note: the gradient has a singularity at x=0 and x=1
            if((x < DataType(0)) || (x > DataType(1)))
              return GradientType::null();
 
            return (_vmax * _den * DataType(4) * (DataType(1) - DataType(2)*x)) * _ve;
          }
 
          HessianType hessian(const PointType& point)
          {
            // project point onto line segment
            const DataType x = Tiny::dot(point - _vo, _ve) * _den;
 
            // Note: the hessian has a singularity at x=0 and x=1
            if((x < DataType(0)) || (x > DataType(1)))
              return HessianType::null();
 
            HessianType hess;
            const DataType v = -DataType(8) * _vmax * _den * _den;
            hess[0][0] = v * _ve[0] * _ve[0];
            hess[0][1] = hess[1][0] = v * _ve[0] * _ve[1];
            hess[1][1] = v * _ve[1] * _ve[1];
            return hess;
          }
        };
      }; // class ParProfileScalar
 
      template<typename DataType_>
      class ParProfileVector :
        public ParProfileBase<DataType_>
      {
      public:
        static constexpr int domain_dim = 2;
        typedef Analytic::Image::Vector<2> ImageType;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
        using ParProfileBase<DataType_>::ParProfileBase;
 
        template<typename Traits_>
        class Evaluator :
          public Analytic::Function::Evaluator<Traits_>
        {
        protected:
          typedef typename Traits_::DataType DataType;
          typedef typename Traits_::PointType PointType;
          typedef typename Traits_::ValueType ValueType;
          typedef typename Traits_::GradientType GradientType;
          typedef typename Traits_::HessianType HessianType;
 
          PointType _vo, _ve, _vn;
          DataType _den, _vmax;
 
        public:
          explicit Evaluator(const ParProfileVector& function)
          {
            _vo[0] = DataType(function._x0);
            _vo[1] = DataType(function._y0);
            _ve[0] = DataType(function._x1 - function._x0);
            _ve[1] = DataType(function._y1 - function._y0);
            _vn[0] =  _ve[1];
            _vn[1] = -_ve[0];
            _vn.normalize();
            _den = DataType(1) / Tiny::dot(_ve, _ve);
            _vmax = DataType(function._vmax);
          }
 
          ValueType value(const PointType& point)
          {
            // project point onto line segment
            const DataType x = Math::clamp(Tiny::dot(point - _vo, _ve) * _den, DataType(0), DataType(1));
            // compute function value
            const DataType v = _vmax * DataType(4) * x * (DataType(1) - x);
            ValueType val;
            val[0] = _vn[0] * v;
            val[1] = _vn[1] * v;
            return val;
          }
 
          GradientType gradient(const PointType& point)
          {
            // project point onto line segment
            const DataType x = Tiny::dot(point - _vo, _ve) * _den;
 
            // Note: the gradient has a singularity at x=0 and x=1
            if((x < DataType(0)) || (x > DataType(1)))
              return GradientType::null();
 
            const DataType v = _vmax * _den * DataType(4) * (DataType(1) - DataType(2)*x);
            GradientType grad;
            grad(0,0) =  v * _vn[0] * _ve[0];
            grad(0,1) =  v * _vn[0] * _ve[1];
            grad(1,0) =  v * _vn[1] * _ve[0];
            grad(1,1) =  v * _vn[1] * _ve[1];
            return grad;
          }
 
          HessianType hessian(const PointType& point)
          {
            // project point onto line segment
            const DataType x = Tiny::dot(point - _vo, _ve) * _den;
 
            // Note: the hessian has a singularity at x=0 and x=1
            if((x < DataType(0)) || (x > DataType(1)))
              return HessianType::null();
 
            const DataType v = -DataType(8) * _vmax * _den * _den;
            HessianType hess;
            hess(0,0,0) = v * _vn[0] * _ve[0] * _ve[0];
            hess(0,1,1) = v * _vn[0] * _ve[1] * _ve[1];
            hess(0,1,0) = hess(0,0,1) = v * _vn[0] * _ve[0] * _ve[1];
            hess(1,1,0) = hess(1,0,1) = v * _vn[1] * _ve[0] * _ve[1];
            hess(1,0,0) = v * _vn[1] * _ve[0] * _ve[0];
            hess(1,1,1) = v * _vn[1] * _ve[1] * _ve[1];
            return hess;
          }
        };
      }; // class ParProfileVector
 
      template<typename DT_>
      class ExpScalarStatic
      {
      public:
        static constexpr DT_ p = DT_(10);
 
        static DT_ eval(DT_ x)
        {
          return (Math::exp(p) - Math::exp(p*x*x)) / (Math::exp(p) - DT_(1));
        }
 
        static DT_ der_x(DT_ x)
        {
          return -DT_(2)*p*x*Math::exp(p*x*x)/(Math::exp(p)-DT_(1));
        }
 
        static DT_ der_xx(DT_ x)
        {
          return -DT_(2)*p*Math::exp(p*x*x)*(DT_(2)*p*x*x+DT_(1))/(Math::exp(p)-DT_(1));
        }
      };
 
      #ifdef FEAT_HAVE_HALFMATH
      template<>
      class ExpScalarStatic<Half>
      {
      public:
        typedef Half DT_;
        static DT_ eval(DT_ x)
        {
          return (Math::exp(DT_(10)) - Math::exp(DT_(10)*x*x)) / (Math::exp(DT_(10)) - DT_(1));
        }
 
        static DT_ der_x(DT_ x)
        {
          return -DT_(2)*DT_(10)*x*Math::exp(DT_(10)*x*x)/(Math::exp(DT_(10))-DT_(1));
        }
 
        static DT_ der_xx(DT_ x)
        {
          return -DT_(2)*DT_(10)*Math::exp(DT_(10)*x*x)*(DT_(2)*DT_(10)*x*x+DT_(1))/(Math::exp(DT_(10))-DT_(1));
        }
      };
      #endif
 
      template<typename DataType_>
      using ExpStatic = Analytic::Common::TensorStatic<ExpScalarStatic<DataType_>, DataType_>;
 
      template<int dim_>
      using ExpFunction = Analytic::StaticWrapperFunction<dim_, ExpStatic, true, true, true>;
 
      template<typename DT_, int dim_>
      class XYPlaneRotation : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<dim_> ImageType;
        static constexpr int domain_dim = dim_;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _angular_velocity;
          const PointType& _origin;
 
        public:
          explicit Evaluator(const XYPlaneRotation& function) :
            _angular_velocity(function._angular_velocity),
            _origin(function._origin)
          {
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val;
            val.format();
            val[0] = _angular_velocity*(-(point[1] - _origin[1]));
            val[1] = _angular_velocity*( (point[0] - _origin[0]));
            return val;
          }
 
          GradientType gradient(const PointType& DOXY(point)) const
          {
            GradientType grad;
            grad.format();
            grad[0][1] = -_angular_velocity;
            grad[1][0] = _angular_velocity;
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
 
        }; // class XYPlaneRotation::Evaluator<...>
 
      public:
        const DataType _angular_velocity;
        const PointType _origin;
 
      public:
        explicit XYPlaneRotation(const DataType angular_velocity, const PointType& origin) :
          _angular_velocity(angular_velocity),
          _origin(origin)
        {
        }
      }; // class XYPlaneRotation
 
      template<typename DT_, int dim_>
      class YZPlaneParabolic : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<dim_> ImageType;
        static constexpr int domain_dim = dim_;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, 2> RangeType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          DataType _fac;
          const std::vector<Tiny::Vector<DataType, 2>>& _range;
 
        public:
          explicit Evaluator(const YZPlaneParabolic& function) :
            _fac(function._amplitude),
            _range(function._range)
 
          {
            XASSERT(_range.size() == size_t(domain_dim-1));
 
            for (int d(1); d < domain_dim; ++d)
            {
              _fac *= DataType(4)/Math::sqr(_range.at(std::size_t(d-1))[1]-_range.at(std::size_t(d-1))[0]);
            }
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val;
            val.format();
            val(0) = _fac;
 
            for (int d(1); d < domain_dim; ++d)
            {
              val(0) *= (point[d] - _range.at(std::size_t(d-1))[0])*(_range.at(std::size_t(d-1))[1] - point[d]);
            }
 
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad;
            grad.format();
            for(int d(1); d < domain_dim; ++d)
            {
              grad[d][0] = _fac*(_range.at(std::size_t(d-1))[0] + _range.at(std::size_t(d-1))[1] - DataType(2)*point(d));
 
              for (int q(1); q < domain_dim; ++q)
              {
                if (q != d)
                  grad[d][0] *= (point[q] - _range.at(std::size_t(q - 1))[0]) * (_range.at(std::size_t(q - 1))[1] - point[q]);
              }
            }
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
 
        }; // class YZPlaneParabolic::Evaluator<...>
 
      public:
        const DataType _amplitude;
        std::vector<Tiny::Vector<DataType, 2>> _range;
 
      public:
        explicit YZPlaneParabolic(const DataType amplitude, const RangeType& range_y) :
          _amplitude(amplitude),
          _range(1)
        {
          _range.at(0) = range_y;
        }
 
        explicit YZPlaneParabolic(const DataType amplitude, const RangeType& range_y, const RangeType& range_z) :
          _amplitude(amplitude),
          _range(2)
        {
          _range.at(0) = range_y;
          _range.at(1) = range_z;
        }
      }; // class YZPlaneParabolic
 
      template<typename DT_, int dim_>
      class SinYT0 : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<dim_> ImageType;
        static constexpr int domain_dim = dim_;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _t;
 
        public:
          explicit Evaluator(const SinYT0& function) :
            _t(function._t)
 
          {
            XASSERT(_t >= DataType(0));
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val;
            val.format();
            val(0) = Math::sin(_t*point[1]);
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad;
            grad.format();
            grad[0][1] = _t*Math::cos(_t*point[1]);
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
 
        }; // class SinYT0::Evaluator<...>
 
      private:
        DataType _t;
 
      public:
        explicit SinYT0(DataType t = DataType(0)) :
          _t(t)
        {
        }
 
        void set_time(const DataType t)
        {
          _t = t;
        }
 
      };
 
      template<typename DT_, int dim_>
      class SinYT0StokesRhs : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<dim_> ImageType;
        static constexpr int domain_dim = dim_;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = false;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _t;
          const DataType _fac;
 
        public:
          explicit Evaluator(const SinYT0StokesRhs& function) :
            _t(function._t),
            _fac(DataType(1)/function._reynolds)
 
          {
            XASSERT(_t >= DataType(0));
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val;
            val.format();
            val(0) = point[1]*Math::cos(_t*point[1]) + _fac*Math::sqr(_t)*Math::sin(point[1]*_t);
            return val;
          }
 
          GradientType gradient(const PointType& DOXY(point)) const
          {
            return GradientType::null();
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
 
        }; // class SinYT0StokesRhs::Evaluator<...>
 
      private:
        const DataType _reynolds;
        DataType _t;
 
      public:
        explicit SinYT0StokesRhs(DataType reynolds, DataType t = DataType(0)) :
          _reynolds(reynolds),
          _t(t)
        {
          XASSERT(reynolds > DataType(0));
        }
 
        void set_time(const DataType t)
        {
          _t = t;
        }
      };
 
      template<typename DT_, int dim_>
      class GuermondStokesSolPressure : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr int domain_dim = dim_;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _t;
 
        public:
          explicit Evaluator(const GuermondStokesSolPressure& function) :
            _t(function._t)
 
          {
            XASSERT(_t >= DataType(0));
          }
 
          ValueType value( const PointType& point)
          {
            ValueType val;
            val = Math::sin(point[0] - point[1] + _t);
            // To make it mean value 0 on [0,1]x[0,1]
            val -= (DataType(2)*Math::sin(_t) - Math::sin(DataType(1)+_t) - Math::sin(-DataType(1) + _t));
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad;
            grad.format();
            grad[0] = Math::cos(point[0] - point[1] + _t);
            grad[1] = -grad[0];
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
        }; // class GuermondStokesSolPressure::Evaluator<...>
 
      private:
        DataType _t;
 
      public:
        explicit GuermondStokesSolPressure(DataType t = DataType(0)) :
          _t(t)
        {
        }
 
        void set_time(const DataType t)
        {
          _t = t;
        }
 
      };
 
      template<typename DT_, int dim_>
      class GuermondStokesSol : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<dim_> ImageType;
        static constexpr int domain_dim = dim_;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _t;
 
        public:
          explicit Evaluator(const GuermondStokesSol& function) :
            _t(function._t)
          {
            XASSERT(_t >= DataType(0));
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val;
            val.format();
            val(0) = Math::sin(point[0] + _t)*Math::sin(point[1] + _t);
            val(1) = Math::cos(point[0] + _t)*Math::cos(point[1] + _t);
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad;
            grad.format();
 
            grad[0][0] = Math::cos(point[0] + _t)*Math::sin(point[1] + _t);
            grad[0][1] = Math::sin(point[0] + _t)*Math::cos(point[1] + _t);
 
            grad[1][0] = -Math::sin(point[0] + _t)*Math::cos(point[1] + _t);
            grad[1][1] = -Math::cos(point[0] + _t)*Math::sin(point[1] + _t);
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
 
        }; // class GuermondStokesSol::Evaluator<...>
 
      private:
        DataType _t;
 
      public:
        explicit GuermondStokesSol(DataType t = DataType(0)) :
          _t(t)
        {
        }
 
        void set_time(const DataType t)
        {
          _t = t;
        }
 
      };
 
      template<typename DT_, int dim_>
      class GuermondStokesSolRhs : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<dim_> ImageType;
        static constexpr int domain_dim = dim_;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = false;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _t;
          const DataType _fac;
 
        public:
          explicit Evaluator(const GuermondStokesSolRhs& function) :
            _t(function._t),
            _fac(DataType(1)/function._reynolds)
 
          {
            XASSERT(_t >= DataType(0));
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val;
            val.format();
 
            // Stationary
            //val(0) = _fac*DataType(2)*Math::sin(point[0]+_t)*Math::sin(point[1]+_t) + Math::cos(point[0]-point[1]+_t);
            //val(1) = _fac*DataType(2)*Math::cos(point[0]+_t)*Math::cos(point[1]+_t) - Math::cos(point[0]-point[1]+_t);
 
            val(0) =  Math::cos(point[0]+_t)*Math::sin(point[1]+_t) + Math::sin(point[0]+_t)*Math::cos(point[1]+_t)
              + _fac*DataType(2)*Math::sin(point[0]+_t)*Math::sin(point[1]+_t) + Math::cos(point[0]-point[1]+_t);
            val(1) = -Math::sin(point[0]+_t)*Math::cos(point[1]+_t) - Math::cos(point[0]+_t)*Math::sin(point[1]+_t)
              + _fac*DataType(2)*Math::cos(point[0]+_t)*Math::cos(point[1]+_t) - Math::cos(point[0]-point[1]+_t);
            return val;
          }
 
          GradientType gradient(const PointType& DOXY(point)) const
          {
            return GradientType::null();
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
 
        }; // class GuermondStokesSolRhs::Evaluator<...>
 
      private:
        const DataType _reynolds;
        DataType _t;
 
      public:
        explicit GuermondStokesSolRhs(DataType reynolds, DataType t = DataType(0)) :
          _reynolds(reynolds),
          _t(t)
        {
          XASSERT(reynolds > DataType(0));
        }
 
        void set_time(const DataType t)
        {
          _t = t;
        }
      };
 
 
      class SphereSinBubbleFunction : public Analytic::Function
      {
      public:
        typedef Analytic::Image::Scalar ImageType;
        static constexpr int domain_dim = 3;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType pi;
 
        public:
          explicit Evaluator(const SphereSinBubbleFunction&) :
            pi(Math::pi<DataType>())
          {
          }
 
          ValueType value(const PointType& point)
          {
            const DataType re = pi / point.norm_euclid();
            return Math::sin(re*point[0]) * Math::sin(re*point[1]) * Math::sin(re*point[2]);
          }
 
          GradientType gradient(const PointType& point) const
          {
            const DataType x = point[0];
            const DataType y = point[1];
            const DataType z = point[2];
            const DataType re = DataType(1) / point.norm_euclid();
            const DataType cx = Math::cos(pi * re * x);
            const DataType cy = Math::cos(pi * re * y);
            const DataType cz = Math::cos(pi * re * z);
            const DataType sx = Math::sin(pi * re * x);
            const DataType sy = Math::sin(pi * re * y);
            const DataType sz = Math::sin(pi * re * z);
            GradientType grad;
            grad[0] = +pi * Math::cub(re) * (cx*sz*sy*y*y + cx*sz*sy*z*z - sx*x*y*cy*sz - sx*sy*x*z*cz);
            grad[1] = -pi * Math::cub(re) * (x*y*cx*sy*sz - sz*sx*cy*x*x - sz*sx*cy*z*z + sx*sy*z*y*cz);
            grad[2] = -pi * Math::cub(re) * (x*z*cx*sy*sz + sx*z*y*cy*sz - sx*cz*sy*x*x - sx*cz*sy*y*y);
            return grad;
          }
 
          HessianType hessian(const PointType& point) const
          {
            const DataType x = point[0];
            const DataType y = point[1];
            const DataType z = point[2];
            const DataType di = Math::sqrt(x*x + y*y + z*z);
            const DataType re = DataType(1) / di;
            const DataType qo = DataType(1) / Math::cub(x*x + y*y + z*z);
            const DataType cx = Math::cos(pi * re * x);
            const DataType cy = Math::cos(pi * re * y);
            const DataType cz = Math::cos(pi * re * z);
            const DataType sx = Math::sin(pi * re * x);
            const DataType sy = Math::sin(pi * re * y);
            const DataType sz = Math::sin(pi * re * z);
            HessianType hess;
            hess[0][0] = -pi*qo*(
              DataType(2)*cx*sy*cz*pi*x*y*y*z +
              DataType(2)*cx*sy*cz*pi*x*z*z*z +
              DataType(2)*cx*sz*cy*pi*x*y*y*y +
              DataType(2)*cx*sz*cy*pi*x*y*z*z +
              DataType(2)*sy*sz*sx*pi*y*y*z*z +
              -DataType(2)*sx*y*x*x*cy*sz*di +
              -DataType(2)*sx*y*pi*x*x*cy*z*cz +
              -DataType(2)*sx*sy*z*x*x*cz*di +
              DataType(3)*cx*sy*sz*x*di*y*y +
              DataType(3)*cx*sy*sz*x*di*z*z +
              sx*y*y*pi*x*x*sy*sz +
              sx*sy*z*z*pi*x*x*sz +
              sy*sz*sx*pi*y*y*y*y +
              sy*sz*sx*pi*z*z*z*z +
              sy*cz*sx*di*y*y*z +
              sy*cz*sx*di*z*z*z +
              sz*cy*sx*di*y*y*y +
              sz*cy*sx*di*y*z*z);
            hess[1][1] = -pi*qo*(
              DataType(2)*cx*sz*cy*pi*x*x*x*y +
              DataType(2)*cx*sz*cy*pi*x*y*z*z +
              DataType(2)*sx*sy*z*z*pi*x*x*sz +
              DataType(2)*sx*y*pi*x*x*cy*z*cz +
              DataType(2)*sx*cy*cz*pi*y*z*z*z +
              -DataType(2)*cx*sy*cz*pi*x*y*y*z +
              -DataType(2)*cx*sy*sz*x*di*y*y +
              -DataType(2)*sy*cz*sx*di*y*y*z +
              DataType(3)*sx*y*x*x*cy*sz*di +
              DataType(3)*sz*cy*sx*di*y*z*z +
              sy*sz*sx*pi*x*x*x*x +
              sx*y*y*pi*x*x*sy*sz +
              sy*sz*sx*pi*y*y*z*z +
              sy*sz*sx*pi*z*z*z*z +
              di*cx*sy*sz*x*x*x +
              cx*sy*sz*x*di*z*z +
              sx*sy*z*x*x*cz*di +
              sy*cz*sx*di*z*z*z);
            hess[2][2] = -pi*qo*(
              DataType(2)*cx*sy*cz*pi*x*x*x*z +
              DataType(2)*cx*sy*cz*pi*x*y*y*z +
              DataType(2)*sx*y*y*pi*x*x*sy*sz +
              DataType(2)*sx*y*pi*x*x*cy*z*cz +
              DataType(2)*sx*cy*cz*pi*y*y*y*z +
              -DataType(2)*cx*sy*sz*x*di*z*z +
              -DataType(2)*cx*sz*cy*pi*x*y*z*z +
              -DataType(2)*sz*cy*sx*di*y*z*z +
              DataType(3)*sx*sy*z*x*x*cz*di +
              DataType(3)*sy*cz*sx*di*y*y*z +
              sy*sz*sx*pi*x*x*x*x +
              sx*sy*z*z*pi*x*x*sz +
              sy*sz*sx*pi*y*y*y*y +
              sy*sz*sx*pi*y*y*z*z +
              di*cx*sy*sz*x*x*x +
              cx*sy*sz*x*di*y*y +
              sx*y*x*x*cy*sz*di +
              sz*cy*sx*di*y*y*y);
            hess[0][1] = hess[1][0] = pi*qo*(
              DataType(2)*y*y*pi*x*x*cx*cy*sz +
              DataType(2)*di*cx*sy*sz*x*x*y +
              DataType(2)*sx*y*y*x*cy*sz*di +
              DataType(3)*sx*sy*z*x*cz*y*di +
              y*pi*x*x*cx*sy*z*cz +
              -cx*sy*cz*pi*y*y*y*z +
              -cx*sy*cz*pi*y*z*z*z +
              cx*sz*cy*pi*x*x*z*z +
              cx*sz*cy*pi*y*y*z*z +
              cx*sz*cy*pi*z*z*z*z +
              sy*sz*sx*pi*x*x*x*y +
              sy*sz*sx*pi*x*y*y*y +
              sx*sy*z*z*pi*x*y*sz +
              -sx*cy*cz*pi*x*x*x*z +
              sx*y*y*pi*x*cy*z*cz +
              -sx*cy*cz*pi*x*z*z*z +
              -cx*sy*sz*di*y*y*y +
              -cx*sy*sz*di*y*z*z +
              -cy*sz*sx*x*x*x*di +
              -cy*sz*sx*x*di*z*z);
            hess[0][2] = hess[2][0] = pi*qo*(
              DataType(2)*z*z*pi*x*x*cx*sy*cz +
              DataType(2)*sx*sy*z*z*x*cz*di +
              DataType(2)*di*cx*sy*sz*x*x*z +
              DataType(3)*sx*y*x*cy*sz*z*di +
              cx*sy*cz*pi*x*x*y*y +
              cx*sy*cz*pi*y*y*y*y +
              cx*sy*cz*pi*y*y*z*z +
              sy*sz*sx*pi*x*x*x*z +
              sx*y*y*pi*x*z*sy*sz +
              sy*sz*sx*pi*x*z*z*z +
              -sx*cy*cz*pi*x*x*x*y +
              -sx*cy*cz*pi*x*y*y*y +
              sx*y*pi*z*z*cy*x*cz +
              z*pi*x*x*cx*y*cy*sz +
              -cx*sz*cy*pi*y*y*y*z +
              -cx*sz*cy*pi*y*z*z*z +
              -sy*cz*sx*x*x*x*di +
              -sy*cz*sx*x*di*y*y +
              -sy*sz*cx*di*y*y*z +
              -sy*sz*cx*di*z*z*z);
            hess[1][2] = hess[2][1] = pi*qo*(
              DataType(2)*sx*y*y*pi*z*z*cy*cz +
              DataType(2)*sx*sy*z*z*y*cz*di +
              DataType(2)*di*sz*sx*cy*y*y*z +
              DataType(3)*y*x*cx*sy*sz*z*di +
              -cx*sy*cz*pi*x*x*x*y +
              -cx*sy*cz*pi*x*y*y*y +
              z*z*pi*x*cx*sy*y*cz +
              -cx*sz*cy*pi*x*x*x*z +
              y*y*pi*x*cx*z*cy*sz +
              -cx*sz*cy*pi*x*z*z*z +
              y*pi*x*x*z*sx*sy*sz +
              sy*sz*sx*pi*y*y*y*z +
              sy*sz*sx*pi*y*z*z*z +
              sx*cy*cz*pi*x*x*x*x +
              sx*cy*cz*pi*x*x*y*y +
              sx*cy*cz*pi*x*x*z*z +
              -sy*cz*sx*x*x*di*y +
              -sy*cz*sx*di*y*y*y +
              -sz*cy*sx*x*x*di*z +
              -sz*cy*sx*di*z*z*z);
            return hess;
          }
        }; // class SphereSinBubbleFunction::Evaluator<...>
      }; // class SphereSinBubbleFunction
 
      class StandingVortexFunction2D :
        public Analytic::Function
      {
      public:
        static constexpr int domain_dim = 2;
        typedef Analytic::Image::Vector<2> ImageType;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = false;
 
        template<typename Traits_>
        class Evaluator :
          public Analytic::Function::Evaluator<Traits_>
        {
        protected:
          typedef typename Traits_::DataType DataType;
          typedef typename Traits_::PointType PointType;
          typedef typename Traits_::ValueType ValueType;
          typedef typename Traits_::GradientType GradientType;
 
        public:
          explicit Evaluator(const StandingVortexFunction2D&)
          {
          }
 
          ValueType value(const PointType& point) const
          {
            ValueType val(DataType(0));
            const DataType x = point[0] - DataType(0.5);
            const DataType y = point[1] - DataType(0.5);
            const DataType r = Math::sqrt(x*x + y*y);
            if (r < DataType(0.2))
            {
              val[0] =  DataType(5)*y;
              val[1] = -DataType(5)*x;
            }
            else if (r < DataType(0.4))
            {
              val[0] = ( DataType(2)/r - DataType(5))*y;
              val[1] = (-DataType(2)/r + DataType(5))*x;
            }
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad(DataType(0));
            const DataType x = point[0] - DataType(0.5);
            const DataType y = point[1] - DataType(0.5);
            const DataType r = Math::sqrt(x*x + y*y);
            if (r < DataType(0.2))
            {
              grad[0][1] =  DataType(5);
              grad[1][0] = -DataType(5);
            }
            else if (r < DataType(0.4))
            {
              const DataType s = DataType(1) / (r*r*r);
              grad[0][0] = -DataType(2)*x*y*s;
              grad[0][1] = -DataType(2)*y*y*s + DataType(2)/r - DataType(5);
              grad[1][0] =  DataType(2)*x*x*s - DataType(2)/r + DataType(5);
              grad[1][1] =  DataType(2)*x*y*s;
            }
            return grad;
          }
        };
      }; // class StandingVortexFunction2D
 
      template<typename DT_>
      class TaylorGreenVortexVelo2D : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<2> ImageType;
        static constexpr int domain_dim = 2;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        DataType nu;
        DataType cur_t;
 
        explicit TaylorGreenVortexVelo2D(DataType nu_ = DataType(1), DataType t_ = DataType(0)) :
          nu(nu_),
          cur_t(t_)
        {
        }
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _pi;
          const DataType _fac;
 
        public:
          explicit Evaluator(const TaylorGreenVortexVelo2D& function) :
            _pi(Math::pi<DataType>()),
            _fac(Math::exp(-DataType(2)*_pi*_pi*function.nu*function.cur_t))
          {
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val = ValueType::null();
            val[0] = +_fac * Math::sin(_pi*point[0]) * Math::cos(_pi*point[1]);
            val[1] = -_fac * Math::cos(_pi*point[0]) * Math::sin(_pi*point[1]);
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad = GradientType::null();
            grad[0][0] = +_fac * _pi * Math::cos(_pi*point[0]) * Math::cos(_pi*point[1]);
            grad[0][1] = -_fac * _pi * Math::sin(_pi*point[0]) * Math::sin(_pi*point[1]);
            grad[1][0] = +_fac * _pi * Math::sin(_pi*point[0]) * Math::sin(_pi*point[1]);
            grad[1][1] = -_fac * _pi * Math::cos(_pi*point[0]) * Math::cos(_pi*point[1]);
            return grad;
          }
 
          HessianType hessian(const PointType& point) const
          {
            HessianType hess = HessianType::null();
            hess[0][0][0] = -_fac * _pi*_pi * Math::sin(_pi*point[0]) * Math::cos(_pi*point[1]);
            hess[0][0][1] = -_fac * _pi*_pi * Math::cos(_pi*point[0]) * Math::sin(_pi*point[1]);
            hess[0][1][0] = -_fac * _pi*_pi * Math::cos(_pi*point[0]) * Math::sin(_pi*point[1]);
            hess[0][1][1] = -_fac * _pi*_pi * Math::sin(_pi*point[0]) * Math::cos(_pi*point[1]);
            hess[1][0][0] = +_fac * _pi*_pi * Math::cos(_pi*point[0]) * Math::sin(_pi*point[1]);
            hess[1][0][1] = +_fac * _pi*_pi * Math::sin(_pi*point[0]) * Math::cos(_pi*point[1]);
            hess[1][1][0] = +_fac * _pi*_pi * Math::sin(_pi*point[0]) * Math::cos(_pi*point[1]);
            hess[1][1][1] = +_fac * _pi*_pi * Math::cos(_pi*point[0]) * Math::sin(_pi*point[1]);
            return hess;
          }
        }; // class TaylorGreenVortexVelo::Evaluator<...>
      }; // class TaylorGreenVortexVelo<...>
 
      template<typename DT_>
      class TaylorGreenVortexPres2D : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr int domain_dim = 2;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        DataType nu;
        DataType cur_t;
 
        explicit TaylorGreenVortexPres2D(DataType nu_ = DataType(1), DataType t_ = DataType(0)) :
          nu(nu_),
          cur_t(t_)
        {
        }
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _pi;
          const DataType _fac;
 
        public:
          explicit Evaluator(const TaylorGreenVortexPres2D& function) :
            _pi(Math::pi<DataType>()),
            _fac(Math::exp(-DataType(4)*_pi*_pi*function.nu*function.cur_t) * DataType(0.5))
          {
          }
 
          ValueType value(const PointType& point)
          {
            return ValueType(_fac * (Math::sqr(Math::cos(_pi*point[0])) + Math::sqr(Math::cos(_pi*point[1])) - DataType(1)));
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad = GradientType::null();
            grad[0] = -DataType(2) * _pi * _fac * Math::sin(_pi*point[0]) * Math::cos(_pi*point[0]);
            grad[1] = -DataType(2) * _pi * _fac * Math::sin(_pi*point[1]) * Math::cos(_pi*point[1]);
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
        }; // class TaylorGreenVortexPres::Evaluator<...>
      }; // class TaylorGreenVortexPres<...>
 
      template<typename DataType_, int dim_>
      class PoiseuillePipeFlow :
        public Analytic::Function
      {
      public:
        typedef DataType_ DataType;
        typedef Analytic::Image::Vector<dim_> ImageType;
        static constexpr int domain_dim = dim_;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
        typedef Tiny::Vector<DataType, domain_dim> PointType;
 
        PointType origin;
        PointType axis;
        DataType radius;
        DataType v_max;
 
        explicit PoiseuillePipeFlow(PointType origin_, PointType axis_, DataType radius_ = DataType(1), DataType v_max_ = DataType(1)) :
          origin(origin_),
          axis(axis_),
          radius(radius_),
          v_max(v_max_)
        {
        }
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const PointType _origin;
          const PointType _axis;
          const DataType _v_max;
          const DataType _radius_s;
 
        public:
          explicit Evaluator(const PoiseuillePipeFlow& function) :
            _origin(function.origin),
            _axis(PointType(function.axis).normalize()),
            _v_max(function.v_max),
            _radius_s(DataType(1) / Math::sqr(function.radius))
          {
          }
 
          ValueType value(const PointType& point)
          {
            // compute ray to project point onto axis
            const PointType q = _origin + Tiny::dot(point - _origin, _axis) * _axis - point;
            // compute velocity based on distance to axis and in direction of axis
            return _v_max * (DataType(1) - q.norm_euclid_sqr()*_radius_s) * _axis;
          }
 
          GradientType gradient(const PointType& point) const
          {
            // compute ray to project point onto axis
            const PointType q = _origin + Tiny::dot(point - _origin, _axis) * _axis - point;
            // compute gradient
            GradientType grad = GradientType::null();
            //grad.add_outer_product(_axis, Tiny::dot(_axis, q)*_axis - q, -DataType(2) * _radius_s * _v_max);
            grad.set_outer_product(_axis, q);
            grad.add_outer_product(_axis, _axis, -Tiny::dot(_axis, q));
            return DataType(2) * _radius_s * _v_max * grad;
          }
 
          HessianType hessian(const PointType&) const
          {
            // seriously, getting this mumbo jumbo below right was mostly a lucky guess...
            GradientType grad_q = GradientType::null();
            grad_q.set_outer_product(_axis, _axis);
            grad_q.add_scalar_main_diag(-DataType(1));
            GradientType at_x_a;
            at_x_a.set_outer_product(_axis, _axis);
            HessianType hess = HessianType::null();
            hess.add_vec_mat_outer_product(_axis, grad_q);
            hess.add_vec_mat_outer_product(_axis * grad_q, at_x_a, -DataType(1));
            return DataType(2) * _radius_s * _v_max * hess;
          }
        }; // class PoiseuillePipeFlow::Evaluator<...>
      }; // class PoiseuillePipeFlow
 
      template<typename DT_>
      class RigidBodyVortexVelo2D : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<2> ImageType;
        static constexpr int domain_dim = 2;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        public:
          explicit Evaluator(const RigidBodyVortexVelo2D&)
          {
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val = ValueType::null();
            val[0] = -point[1];
            val[1] =  point[0];
            return val;
          }
 
          GradientType gradient(const PointType& DOXY(point)) const
          {
            GradientType grad = GradientType::null();
            grad[0][1] = -DataType(1);
            grad[1][0] =  DataType(1);
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
        }; // class RigidBodyVortexVelo2D::Evaluator<...>
      }; // class RigidBodyVortexVelo2D<...>
 
      template<typename DT_>
      class RigidBodyVortexPres2D : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr int domain_dim = 2;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        public:
          explicit Evaluator(const RigidBodyVortexPres2D&)
          {
          }
 
          ValueType value(const PointType& point)
          {
            return DataType(0.5) * (point[0]*point[0] + point[1]*point[1]) - DataType(0.25);
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad = GradientType::null();
            grad[0] = point[0];
            grad[1] = point[1];
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            HessianType hess = HessianType::null();
            hess[0][0] = DataType(1);
            hess[1][1] = DataType(1);
            return hess;
          }
        }; // class RigidBodyVortexPres2D::Evaluator<...>
      }; // class RigidBodyVortexPres2D<...>
 
      template<typename DT_>
      class SineRingVortexVelo2D : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<2> ImageType;
        static constexpr int domain_dim = 2;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _pi2;
 
        public:
          explicit Evaluator(const SineRingVortexVelo2D&) :
            _pi2(DataType(2) * Math::pi<DataType>())
          {
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val = ValueType::null();
            const DataType d = point.norm_euclid();
            if(d < DataType(1E-3))
              return val;
 
            const DataType s2pdi = Math::sin(_pi2 * d) / d;
 
            val[0] =  point[1]*s2pdi;
            val[1] = -point[0]*s2pdi;
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad = GradientType::null();
            const DataType d = point.norm_euclid();
            if(d < DataType(1E-3))
              return grad;
 
            const DataType x(point[0]), y(point[1]);
            const DataType di2 = DataType(1) / (d*d);
            const DataType s2pdi = Math::sin(_pi2 * d) / d;
            const DataType c2pdi2 = Math::cos(_pi2 * d) * di2;
 
            grad[0][0] = -y*x*s2pdi*di2 + _pi2*y*x*c2pdi2;
            grad[0][1] = -y*y*s2pdi*di2 + _pi2*y*y*c2pdi2 + s2pdi;
            grad[1][0] =  x*x*s2pdi*di2 - _pi2*x*x*c2pdi2 - s2pdi;
            grad[1][1] =  x*y*s2pdi*di2 - _pi2*y*x*c2pdi2;
            return grad;
          }
 
          HessianType hessian(const PointType& point) const
          {
            HessianType hess = HessianType::null();
            const DataType d = point.norm_euclid();
            if(d < DataType(1E-3))
              return hess;
 
            const DataType x(point[0]), y(point[1]);
            const DataType di2 = DataType(1) / (d*d);
            const DataType s2pdi3 = Math::sin(_pi2 * d) * di2 / d;
            const DataType c2pdi2 = Math::cos(_pi2 * d) * di2;
            const DataType t3 = DataType(3);
 
            hess[0][0][0] =  t3*y*x*x*s2pdi3*di2 - t3*_pi2*y*x*x*c2pdi2*di2 - y*s2pdi3 + _pi2*y*c2pdi2 - _pi2*_pi2*y*x*x*s2pdi3;
            hess[0][0][1] =
            hess[0][1][0] =  t3*y*y*x*s2pdi3*di2 - x*s2pdi3 - t3*_pi2*y*y*x*c2pdi2*di2 + _pi2*x*c2pdi2 - _pi2*_pi2*y*y*x*s2pdi3;
            hess[0][1][1] =  t3*y*y*y*s2pdi3*di2 - t3*y*s2pdi3 - t3*_pi2*y*y*y*c2pdi2*di2 + t3*_pi2*y*c2pdi2 - _pi2*_pi2*y*y*y*s2pdi3;
            hess[1][0][0] = -t3*x*x*x*s2pdi3*di2 + t3*x*s2pdi3 + t3*_pi2*x*x*x*c2pdi2*di2 - t3*_pi2*x*c2pdi2 + _pi2*_pi2*x*x*x*s2pdi3;
            hess[1][0][1] =
            hess[1][1][0] = -t3*y*x*x*s2pdi3*di2 + t3*_pi2*y*x*x*c2pdi2*di2 + y*s2pdi3 - _pi2*y*c2pdi2 + _pi2*_pi2*y*x*x*s2pdi3;
            hess[1][1][1] = -t3*y*y*x*s2pdi3*di2 + x*s2pdi3 + t3*_pi2*y*y*x*c2pdi2*di2 - _pi2*x*c2pdi2 + _pi2*_pi2*y*y*x*s2pdi3;
            return hess;
          }
        }; // class SineRingVortexVelo2D::Evaluator<...>
      }; // class SineRingVortexVelo2D<...>
 
      template<typename DT_>
      class SineRingVortexPres2D : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr int domain_dim = 2;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        DataType mean_shift = DataType(1) / Math::sqr(Math::pi<DataType>());
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _pi2, _spi;
 
        public:
          explicit Evaluator(const SineRingVortexPres2D& func) :
            _pi2(DataType(2) * Math::pi<DataType>()),
            _spi(func.mean_shift)
          {
          }
 
          ValueType value(const PointType& point)
          {
            const DataType d = point.norm_euclid();
            return -Math::sin(_pi2*d) / _pi2 - _spi;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad = GradientType::null();
            const DataType d = point.norm_euclid();
            if(d < DataType(1E-3))
              return grad;
 
            const DataType c2pdi = Math::cos(_pi2 * d)  / d;
 
            grad[0] = -point[0]*c2pdi;
            grad[1] = -point[1]*c2pdi;
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
        }; // class SineRingVortexPres2D::Evaluator<...>
      }; // class SineRingVortexPres2D<...>
 
      template<typename DT_>
      class SineRingVortexRHS2D : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Vector<2> ImageType;
        static constexpr int domain_dim = 2;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = false;
        static constexpr bool can_hess = false;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        DataType nu, beta, theta, sigma;
 
        explicit SineRingVortexRHS2D(DataType nu_ = DataType(1), DataType beta_ = DataType(0), DataType theta_ = DataType(0), DataType sigma_ = DataType(1)) :
          nu(nu_),
          beta(beta_),
          theta(theta_),
          sigma(sigma_)
        {
        }
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const DataType _pi2, _nu, _beta, _theta, _sigma;
 
        public:
          explicit Evaluator(const SineRingVortexRHS2D& function) :
            _pi2(DataType(2) * Math::pi<DataType>()),
            _nu(function.nu),
            _beta(function.beta),
            _theta(function.theta),
            _sigma(function.sigma)
          {
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val = ValueType::null();
            const DataType d = point.norm_euclid();
            if(d < DataType(1E-3))
              return val;
 
            const DataType x(point[0]), y(point[1]);
            const DataType q = DataType(1) / d;
            const DataType s2pdq = Math::sin(_pi2 * d) * q;
            const DataType c2pdq = Math::cos(_pi2 * d) * q;
 
            val[0] = -_nu*q*y*( _pi2*c2pdq - (_pi2*_pi2*d + q)*s2pdq) - _beta*x*s2pdq*s2pdq + _theta*y*s2pdq - _sigma*x*c2pdq;
            val[1] = -_nu*q*x*(-_pi2*c2pdq + (_pi2*_pi2*d + q)*s2pdq) - _beta*y*s2pdq*s2pdq - _theta*x*s2pdq - _sigma*y*c2pdq;
            return val;
          }
        }; // class SineRingVortexRHS2D::Evaluator<...>
      }; // class SineRingVortexRHS2D<...>
 
      template<int dim_, int a_ = 2>
      class SphereCapFunction :
        public Analytic::Function
      {
      public:
        static_assert(a_ >= 1, "parameter a_ must be >= 1");
 
        static constexpr int domain_dim = dim_;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
        template<typename Traits_>
        class Evaluator :
          public Analytic::Function::Evaluator<Traits_>
        {
        public:
          typedef typename Traits_::DataType DataType;
          typedef typename Traits_::PointType PointType;
          typedef typename Traits_::ValueType ValueType;
          typedef typename Traits_::GradientType GradientType;
          typedef typename Traits_::HessianType HessianType;
 
          const DataType sa;
 
          explicit Evaluator(const SphereCapFunction&) :
            sa(Math::sqrt(DataType(a_) - DataType(1)))
          {
          }
 
          ValueType value(const PointType& point) const
          {
            return Math::sqrt(DataType(a_) - point.norm_euclid_sqr()) - sa;
          }
 
          GradientType gradient(const PointType& point) const
          {
            return (DataType(-1) / Math::sqrt(DataType(a_) - point.norm_euclid_sqr())) * point;
          }
 
          HessianType hessian(const PointType& point) const
          {
            const DataType s = DataType(1) / Math::cub(Math::sqrt(DataType(a_) - point.norm_euclid_sqr()));
            //const DataType s = DataType(1) / Math::pow(a - point.norm_euclid_sqr(), DataType(1.5));
            HessianType hess;
            hess.format();
            for(int i = 0; i < dim_; ++i)
            {
              for(int j = 0; j < dim_; ++j)
              {
                if(i == j)
                {
                  hess[i][j] = DataType(0);
                  for(int k = 0; k < dim_; ++k)
                    hess[i][j] += (k == i ? DataType(-a_) : Math::sqr(point[k])) * s;
                }
                else
                {
                  hess[i][j] = -point[i] * point[j] * s;
                }
              }
            }
            return hess;
          }
        }; // class SphereCapFunction::Evaluator<...>
      }; // class SphereCapFunction
 
      template<typename DT_>
      class CornerSingularity2D :
        public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        static constexpr int domain_dim = 2;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
      private:
        Tiny::Vector<DataType, 2> base_point;
        DataType alpha, k, offset;
 
      public:
        CornerSingularity2D() = delete;
 
        CornerSingularity2D(const Tiny::Vector<DT_, 2>& center, Tiny::Vector<DT_, 2> v_1, Tiny::Vector<DT_, 2> v_2, DT_ alpha_ = DT_(1.)) :
          base_point(center),
          alpha(alpha_)
        {
          ASSERTM((v_1.norm_euclid_sqr() > Math::eps<DataType>()) && (v_2.norm_euclid_sqr() > Math::eps<DataType>()), "Spanvectors are zero length");
          k = Math::pi<DataType>() / Tiny::calculate_opening_angle(v_1,v_2);
          offset = Tiny::calculate_opening_angle(Tiny::Vector<DT_, 2>{DataType(1.0), DataType(0.)}, v_1);
        }
 
        template<typename Traits_>
        class Evaluator :
          public Analytic::Function::Evaluator<Traits_>
        {
        public:
          typedef typename Traits_::DataType DataType;
          typedef typename Traits_::PointType PointType;
          typedef typename Traits_::ValueType ValueType;
          typedef typename Traits_::GradientType GradientType;
          typedef typename Traits_::HessianType HessianType;
 
 
          Tiny::Vector<DataType, 2> base_point;
          const DataType alpha, k, offset;
          const DataType eps;
 
          explicit Evaluator(const CornerSingularity2D& func) :
            base_point(func.base_point),
            alpha(func.alpha),
            k(func.k),
            offset(func.offset),
            eps(Math::eps<DataType>())
          {
          }
 
          ValueType value(const PointType& point) const
          {
            //first calculate radius
            const DataType x = point[0] - base_point[0];
            const DataType y = point[1] - base_point[1];
            const DataType r = Math::sqrt(Math::sqr(x) + Math::sqr(y));
            //TODO: Set x_rel to zero if r is very small?? <- Better conditioning?
            const DataType x_rel = r > DataType(10.)* eps ? x/r : DataType(0);
            DataType phi = Math::acos(x_rel);
            phi = (y >= DataType(0)) ? phi : DataType(2)*Math::pi<DataType>()-phi;
            return alpha * Math::pow(r, k) * Math::sin(k * (phi - offset));
          }
 
          GradientType gradient(const PointType& point) const
          {
            const DataType x = point[0] - base_point[0];
            const DataType y = point[1] - base_point[1];
            const DataType r = Math::sqrt(Math::sqr(x) + Math::sqr(y));
            const DataType x_rel = r > DataType(10.)* eps ? x/r : DataType(0);
            const DataType y_rel = r > DataType(10.)* eps ? y/r : DataType(0);
            //Note: While in practice we use a rotated coordinate system to express our polar-coordinate
            //based function, we have to use the "real" cos(phi) and sin(phi) values
            //for our base transformation, which are conveniently x_rel and y_rel
            //TODO: This could be badly conditioned for r -> 0, have to test this...
            DataType phi = Math::acos(x_rel);
            phi = (y >= DataType(0)) ? phi : DataType(2)*Math::pi<DataType>()-phi;
            const DataType ex = k * Math::pow(r, k-DataType(1));
            const DataType sval = Math::sin(k * (phi - offset));
            const DataType cval = Math::cos(k * (phi - offset));
            GradientType grad;
            grad[0] = alpha * ex * (x_rel * sval - y_rel * cval);
            grad[1] = alpha * ex * (y_rel * sval + x_rel * cval);
            return grad;
          }
 
          HessianType hessian(const PointType& point) const
          {
            const DataType x = point[0] - base_point[0];
            const DataType y = point[1] - base_point[1];
            const DataType r = Math::sqrt(Math::sqr(x) + Math::sqr(y));
            const DataType x_rel = r > DataType(10.)* eps ? x/r : DataType(0);
            const DataType y_rel = r > DataType(10.)* eps ? y/r : DataType(0);
            DataType phi = Math::acos(x_rel);
            phi = (y >= DataType(0)) ? phi : DataType(2)*Math::pi<DataType>()-phi;
            const DataType ex = k * (k - DataType(1)) * Math::pow(r, k-DataType(2));
            const DataType sval = Math::sin(k * (phi - offset));
            const DataType cval = Math::cos(k * (phi - offset));
            HessianType hess;
            hess[0][0] = alpha * ex * ((Math::sqr(x_rel) - Math::sqr(y_rel)) * sval - DataType(2) * x_rel * y_rel * cval);
            hess[1][1] = -hess[0][0];
            hess[0][1] = hess[1][0] = alpha * ex * (DataType(2) * y_rel * x_rel * sval + (Math::sqr(x_rel) - Math::sqr(y_rel)) * cval);
            return hess;
          }
        }; // class CornerSingularity2D::Evaluator<...>
      }; // CornerSingularity2D
 
      template<typename DT_>
      class CornerSingularity2DRadial :
        public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        static constexpr int domain_dim = 2;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
      private:
        DataType alpha, k;
 
      public:
        CornerSingularity2DRadial() = delete;
 
        CornerSingularity2DRadial(DT_ theta_, DT_ alpha_ = DT_(1.)) :
          alpha(alpha_),
          k(Math::pi<DataType>() / theta_)
        {}
 
        CornerSingularity2DRadial(const Tiny::Vector<DT_, 2>& vec_r, const Tiny::Vector<DT_, 2>& vec_l, DT_ alpha_ = DT_(1.)) :
          alpha(alpha_),
          k(Math::pi<DataType>() / Math::calc_opening_angle(vec_r[0], vec_r[1], vec_l[0], vec_l[1]))
        {}
 
        template<typename Traits_>
        class Evaluator :
          public Analytic::Function::Evaluator<Traits_>
        {
        public:
          typedef typename Traits_::DataType DataType;
          typedef typename Traits_::PointType PointType;
          typedef typename Traits_::ValueType ValueType;
          typedef typename Traits_::GradientType GradientType;
          typedef typename Traits_::HessianType HessianType;
 
 
          const DataType alpha, k;
 
          explicit Evaluator(const CornerSingularity2DRadial& func) :
            alpha(func.alpha),
            k(func.k)
          {
          }
 
          ValueType value(const PointType& point) const
          {
            return alpha * Math::pow(point[0], k) * Math::sin(k * point[1]);
          }
 
          GradientType gradient(const PointType& point) const
          {
            const DataType ex = k * Math::pow(point[0], k-DataType(1));
            const DataType sval = Math::sin(k * point[1]);
            const DataType cval = Math::cos(k * point[1]);
            GradientType grad;
            grad[0] = alpha * ex * sval;
            grad[1] = alpha * ex * cval * point[0];
            return grad;
          }
 
          HessianType hessian(const PointType& point) const
          {
            const DataType ex = k * Math::pow(point[0], k-DataType(2));
            const DataType sval = Math::sin(k * point[1]);
            const DataType cval = Math::cos(k * point[1]);
            HessianType hess;
            hess[0][0] = alpha * ex * (k-DataType(1)) * sval;
            hess[1][1] = - alpha * ex * Math::sqr(point[0]) * k * sval;
            hess[0][1] = hess[1][0] = alpha * ex * k * point[0] * cval;
            return hess;
          }
        }; // class CornerSingularity2DRadial::Evaluator<...>
 
      }; // class CornerSingularity2DRadial
 
      template<typename DT_>
      using CornerSingluratity2DSimple = PolarCoordinate<CornerSingularity2DRadial<DT_>>;
 
      template<typename DT_>
      class FrankesFunction : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef typename Analytic::Image::Scalar ImageType;
        static constexpr int domain_dim = 2;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        public:
          explicit Evaluator(const FrankesFunction&)
          {
          }
 
          ValueType value(const PointType& point)
          {
            const DataType x = point[0];
            const DataType y = point[1];
 
            DataType val = (DataType(3) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)))) / DataType(4) - Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) / DataType(5) +(DataType(3) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) -
                  Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)))) / DataType(4) + Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) / DataType(2);
 
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad;
            const DataType x = point[0];
            const DataType y = point[1];
            grad.format();
            grad[0] = (Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) *(DataType(162) * x - DataType(72))) / DataType(5) -(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(9))) / DataType(4) -
                      (Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(63) / DataType(2))) / DataType(2) -(DataType(3) * Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) /
                      DataType(10)) *((DataType(162) * x) / DataType(49) + DataType(18) / DataType(49))) / DataType(4);
            grad[1] = (Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) *(DataType(162) * y - DataType(126))) / DataType(5) -(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)) *((DataType(81) * y) / DataType(2) - DataType(9))) / DataType(4) -
                      (Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) *((DataType(81) * y) / DataType(2) - DataType(27) / DataType(2))) / DataType(2) -(DataType(27) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) /
                      DataType(10)))) / DataType(40);
            return grad;
          }
 
          HessianType hessian(const PointType& point) const
          {
            HessianType hess;
            const DataType x = point[0];
            const DataType y = point[1];
            hess.format();
            hess[0][0] = (DataType(162) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))))) / DataType(5) -(DataType(243) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)))) / DataType(98) -(DataType(243) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) /
                          DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)))) / DataType(8) -(DataType(81) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)))) / DataType(4) +(DataType(3) * Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) /
                          DataType(49) - DataType(1) / DataType(10)) * Math::sqr((DataType(162) * x) / DataType(49) + DataType(18) / DataType(49))) / DataType(4) +(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)) * Math::sqr((DataType(81) * x) / DataType(2) - DataType(9))) / DataType(4) +(Math::exp(-
                          Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) * Math::sqr((DataType(81) * x) / DataType(2) - DataType(63) / DataType(2))) / DataType(2) -(Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) * Math::sqr(DataType(162) * x - DataType(72))) / DataType(5);
            hess[0][1] = hess[1][0] = (DataType(27) * Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)) *((DataType(162) * x) / DataType(49) + DataType(18) / DataType(49))) / DataType(40) +(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) /
                          DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(9)) *((DataType(81) * y) / DataType(2) - DataType(9))) / DataType(4) +(Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(63) / DataType(2)) *((DataType(81) *
                          y) / DataType(2) - DataType(27) / DataType(2))) / DataType(2) -(Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) *(DataType(162) * x - DataType(72)) *(DataType(162) * y - DataType(126))) / DataType(5);
            hess[1][1] = (DataType(243) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)))) / DataType(400) +(DataType(162) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))))) / DataType(5) -(DataType(243) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) /
                          DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)))) / DataType(8) -(DataType(81) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)))) / DataType(4) +(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y -
                          DataType(2)) / DataType(4)) * Math::sqr((DataType(81) * y) / DataType(2) - DataType(9))) / DataType(4) +(Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) * Math::sqr((DataType(81) * y) / DataType(2) - DataType(27) / DataType(2))) / DataType(2) -(Math::exp(- Math::sqr(DataType(9) * x -
                          DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) * Math::sqr(DataType(162) * y - DataType(126))) / DataType(5);
            return hess;
          }
 
        }; // class FrankesFunction::Evaluator<...>
 
      public:
        explicit FrankesFunction() = default;
      }; // class FrankesFunction
 
      template<typename DT_>
      class Frankes3DVariantFunction : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef typename Analytic::Image::Scalar ImageType;
        static constexpr int domain_dim = 3;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
        typedef Tiny::Vector<DT_, domain_dim> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        public:
          explicit Evaluator(const Frankes3DVariantFunction&)
          {
          }
 
          ValueType value(const PointType& point)
          {
            const DataType x = point[0];
            const DataType y = point[1];
            const DataType z = point[2];
 
            DataType val = (DataType(3) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)))) / DataType(4) - Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) / DataType(5) +(DataType(3) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) -
                  Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)))) / DataType(4) + Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) / DataType(2);
            val *= Math::exp(-Math::sqr(DataType(2)*z-DataType(0.4)));
 
            return val;
          }
 
          GradientType gradient(const PointType& point) const
          {
            GradientType grad;
            const DataType x = point[0];
            const DataType y = point[1];
            const DataType z = point[2];
 
            const DataType franks_val = (DataType(3) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)))) / DataType(4) - Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) / DataType(5) +(DataType(3) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) -
                  Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)))) / DataType(4) + Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) / DataType(2);
 
 
            const DataType zexp = Math::exp(-Math::sqr(DataType(2)*z-DataType(0.4)));
            const DataType dzexp = -DataType(4)*(DataType(2)*z - DataType(0.4)) * zexp;
 
            grad.format();
            grad[0] = zexp * ((Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) *(DataType(162) * x - DataType(72))) / DataType(5) -(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(9))) / DataType(4) -
                      (Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(63) / DataType(2))) / DataType(2) -(DataType(3) * Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) /
                      DataType(10)) *((DataType(162) * x) / DataType(49) + DataType(18) / DataType(49))) / DataType(4));
            grad[1] = zexp * ((Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) *(DataType(162) * y - DataType(126))) / DataType(5) -(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)) *((DataType(81) * y) / DataType(2) - DataType(9))) / DataType(4) -
                      (Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) *((DataType(81) * y) / DataType(2) - DataType(27) / DataType(2))) / DataType(2) -(DataType(27) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) /
                      DataType(10)))) / DataType(40));
            grad[2] = dzexp * franks_val;
            return grad;
          }
 
          HessianType hessian(const PointType& point) const
          {
            HessianType hess;
            const DataType x = point[0];
            const DataType y = point[1];
            const DataType z = point[2];
            const DataType zexp = Math::exp(-Math::sqr(DataType(2)*z-DataType(0.4)));
            const DataType dzexp = -DataType(4)*(DataType(2)*z - DataType(0.4)) * zexp;
            const DataType franks_val = (DataType(3) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)))) / DataType(4) - Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) / DataType(5) +(DataType(3) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) -
                  Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)))) / DataType(4) + Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) / DataType(2);
            hess.format();
            hess[0][0] = zexp * ((DataType(162) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))))) / DataType(5) -(DataType(243) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)))) / DataType(98) -(DataType(243) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) /
                          DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)))) / DataType(8) -(DataType(81) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)))) / DataType(4) +(DataType(3) * Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) /
                          DataType(49) - DataType(1) / DataType(10)) * Math::sqr((DataType(162) * x) / DataType(49) + DataType(18) / DataType(49))) / DataType(4) +(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)) * Math::sqr((DataType(81) * x) / DataType(2) - DataType(9))) / DataType(4) +(Math::exp(-
                          Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) * Math::sqr((DataType(81) * x) / DataType(2) - DataType(63) / DataType(2))) / DataType(2) -(Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) * Math::sqr(DataType(162) * x - DataType(72))) / DataType(5));
            hess[0][1] = hess[1][0] = zexp * ((DataType(27) * Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)) *((DataType(162) * x) / DataType(49) + DataType(18) / DataType(49))) / DataType(40) +(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) /
                          DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(9)) *((DataType(81) * y) / DataType(2) - DataType(9))) / DataType(4) +(Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(63) / DataType(2)) *((DataType(81) *
                          y) / DataType(2) - DataType(27) / DataType(2))) / DataType(2) -(Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) *(DataType(162) * x - DataType(72)) *(DataType(162) * y - DataType(126))) / DataType(5));
            hess[1][1] = zexp * ((DataType(243) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) / DataType(10)))) / DataType(400) +(DataType(162) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))))) / DataType(5) -(DataType(243) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) /
                          DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)))) / DataType(8) -(DataType(81) *(Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)))) / DataType(4) +(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y -
                          DataType(2)) / DataType(4)) * Math::sqr((DataType(81) * y) / DataType(2) - DataType(9))) / DataType(4) +(Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) * Math::sqr((DataType(81) * y) / DataType(2) - DataType(27) / DataType(2))) / DataType(2) -(Math::exp(- Math::sqr(DataType(9) * x -
                          DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) * Math::sqr(DataType(162) * y - DataType(126))) / DataType(5));
            hess[0][2] = hess[2][0] = dzexp * ((Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) *(DataType(162) * x - DataType(72))) / DataType(5) -(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(9))) / DataType(4) -
                          (Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) *((DataType(81) * x) / DataType(2) - DataType(63) / DataType(2))) / DataType(2) -(DataType(3) * Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) /
                          DataType(10)) *((DataType(162) * x) / DataType(49) + DataType(18) / DataType(49))) / DataType(4));
            hess[1][2] = hess[2][1] = dzexp * ((Math::exp(- Math::sqr(DataType(9) * x - DataType(4)) - Math::sqr(DataType(9) * y - DataType(7))) *(DataType(162) * y - DataType(126))) / DataType(5) -(DataType(3) * Math::exp(- Math::sqr(DataType(9) * x - DataType(2)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(2)) / DataType(4)) *((DataType(81) * y) / DataType(2) - DataType(9))) / DataType(4) -
                      (Math::exp(- Math::sqr(DataType(9) * x - DataType(7)) / DataType(4) - Math::sqr(DataType(9) * y - DataType(3)) / DataType(4)) *((DataType(81) * y) / DataType(2) - DataType(27) / DataType(2))) / DataType(2) -(DataType(27) *(Math::exp(-(DataType(9) * y) / DataType(10) - Math::sqr(DataType(9) * x + DataType(1)) / DataType(49) - DataType(1) /
                      DataType(10)))) / DataType(40));
            hess[2][2] = franks_val * (Math::sqr(DataType(4)*(DataType(2)*z - DataType(0.4))) - DataType(8)) * zexp;
            return hess;
          }
 
        }; // class Frankes3DVariantFunction::Evaluator<...>
 
      public:
        explicit Frankes3DVariantFunction() = default;
      }; // class FrankesFunction
 
      template<int dim_, typename DataType_>
      class HarmonicShellFunction :
        public Analytic::Function
      {
      public:
        static_assert((1 <= dim_) && (dim_ <= 3), "harmonic shell function is only implemented in 1D, 2D and 3D");
 
        static constexpr int domain_dim = dim_;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr bool can_value = true;
        static constexpr bool can_grad = true;
        static constexpr bool can_hess = true;
 
        typedef Tiny::Vector<DataType_, dim_> PointType;
 
        template<typename EvalTraits_>
        class Evaluator :
          public Analytic::Function::Evaluator<EvalTraits_>
        {
        public:
          typedef typename EvalTraits_::DataType DataType;
          typedef typename EvalTraits_::PointType PointType;
          typedef typename EvalTraits_::ValueType ValueType;
          typedef typename EvalTraits_::GradientType GradientType;
          typedef typename EvalTraits_::HessianType HessianType;
 
        private:
          const PointType _origin;
          const DataType _eps;
          const DataType _ln2; // = ln(2)
          const DataType _i2ln2; // = 1 /(2*ln(2))
 
        public:
          explicit Evaluator(const HarmonicShellFunction& function) :
            _origin(function._origin),
            _eps(Math::eps<DataType>()),
            _ln2(Math::log(DataType(2))),
            _i2ln2(DataType(0.5)/_ln2)
          {
          }
 
          ValueType value(const PointType& point) const
          {
            if constexpr (dim_ == 1)
              return DataType(3)  - DataType(2) * Math::abs(point[0] - _origin[0]);
            else if constexpr(dim_ == 2)
              return DataType(1) - _i2ln2 * Math::log((point - _origin).norm_euclid_sqr());
            else if constexpr(dim_ == 3)
              return DataType(1) / (point - _origin).norm_euclid();
            else
              return DataType(0);
          }
 
          GradientType gradient(const PointType& point) const
          {
            const DataType norm_sqr = (point - _origin).norm_euclid_sqr();
            if(norm_sqr < _eps)
              return GradientType::null();
 
            GradientType grad(DataType(0));
            if constexpr (dim_ == 1)
              grad[0] = DataType(-2) * Math::signum(point[0] - _origin[0]);
            else if constexpr(dim_ == 2)
              grad = (DataType(-1) / (_ln2 * norm_sqr)) * (point - _origin);
            else if constexpr(dim_ == 3)
              grad = (DataType(1) / Math::pow(norm_sqr, DataType(1.5))) * (_origin - point);
 
            return grad;
          }
 
          HessianType hessian(const PointType& point) const
          {
            const DataType norm_sqr = (point - _origin).norm_euclid_sqr();
            if(norm_sqr < _eps)
              return HessianType::null();
 
            HessianType hess(DataType(0));
            if constexpr(dim_ == 1)
            {
              // 1D hessian is null
            }
            else if constexpr(dim_ == 2)
            {
              const DataType denom = DataType(1) / (_ln2 * Math::sqr(norm_sqr));
              hess[0][0] = denom * (Math::sqr(point[0] - _origin[0]) - Math::sqr(point[1] - _origin[1]));
              hess[1][1] = denom * (Math::sqr(point[1] - _origin[1]) - Math::sqr(point[0] - _origin[0]));
              hess[0][1] = hess[1][0] = denom * DataType(2) * (point[0] - _origin[0]) * (point[1] - _origin[1]);
            }
            else if constexpr(dim_ == 3)
            {
              const DataType denom = DataType(1) / Math::pow(norm_sqr, DataType(2.5));
              for(int i(0); i < dim_; ++i)
              {
                for(int j(0); j < dim_; ++j)
                {
                  if(i == j)
                  {
                    for(int k(0); k < dim_; ++k)
                      hess[i][j] += DataType(k == i ? 2 : -1) * Math::sqr(point[k] - _origin[k]) * denom;
                  }
                  else
                  {
                    hess[i][j] = DataType(3) * (point[i] - _origin[i]) * (point[j] - _origin[j]) * denom;
                  }
                }
              }
            }
            return hess;
          }
        }; // class HarmonicShellFunction::Evaluator<...>
 
      public:
        PointType _origin;
 
      public:
        HarmonicShellFunction()
        {
          _origin.format();
        }
 
        explicit HarmonicShellFunction(const PointType& origin) :
          _origin(origin)
        {
        }
 
        void set_point(const PointType& origin)
        {
          _origin = origin;
        }
      }; // class HarmonicShellFunction
 
      template<int dim_, typename DataType_>
      class DFG95SteadyInflowFunction:
        public Analytic::Function
      {
      public:
        static constexpr int domain_dim = dim_;
        typedef Analytic::Image::Vector<dim_> ImageType;
        static constexpr bool can_value = true;
 
        DataType_ v_max;
 
      public:
        explicit DFG95SteadyInflowFunction(DataType_ vmax) :
          v_max(vmax)
        {
        }
 
        template<typename Traits_>
        class Evaluator :
          public Analytic::Function::Evaluator<Traits_>
        {
        protected:
          typedef typename Traits_::DataType DataType;
          typedef typename Traits_::PointType PointType;
          typedef typename Traits_::ValueType ValueType;
 
          const DataType _d, _scal;
 
        public:
          explicit Evaluator(const DFG95SteadyInflowFunction& function) :
            _d(DataType(0.41)),
            _scal(DataType(function.v_max) / Math::pow(_d, DataType_(2*(dim_-1))))
          {
          }
 
          ValueType value(const PointType& point)
          {
            ValueType val;
            val[0] = _scal;
            for(int i = 1; i < dim_; ++i)
            {
              val[0] *= DataType(4) * point[i] * (_d - point[i]);
              val[i] = DataType(0);
            }
            return val;
          }
        };
      }; // class DFG95SteadyInflowFunction
    } // namespace Common
  } // namespace Analytic
} // namespace FEAT