distance_function.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>
 
namespace FEAT
{
  namespace Analytic
  {
    namespace Distance
    {
      template<int dim_, typename DataType_>
      class DistanceFunction :
        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;
 
        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;
 
        public:
          explicit Evaluator(const DistanceFunction& function) :
            _origin(function._origin)
          {
          }
 
          ValueType value(const PointType& point) const
          {
            return (point - _origin).norm_euclid();
          }
 
          GradientType gradient(const PointType& point) const
          {
            DataType norm (this->value(point));
 
            if(norm < Math::eps<DataType>())
              return GradientType::null();
 
            return (DataType(1) / norm) * (point - _origin);
          }
 
          HessianType hessian(const PointType& point) const
          {
            DataType norm (this->value(point));
            if(norm < Math::eps<DataType>())
              return HessianType::null();
 
            HessianType hess (DataType(0));
            norm = DataType(1)/norm;
            DataType denom = Math::sqr(norm)*norm;
 
            for(int i(0); i < dim_; ++i)
            {
              hess[i][i] = norm;
              for(int j(0); j < dim_; ++j)
              {
                hess[i][j] -= ((point[i] - _origin[i]) * (point[j] - _origin[j]))*denom;
              }
            }
            return hess;
          }
        }; // class DistanceFunction::Evaluator<...>
 
      public:
        PointType _origin;
 
      public:
        DistanceFunction()
        {
          _origin.format();
        }
 
        explicit DistanceFunction(const PointType& origin) :
          _origin(origin)
        {
        }
 
        void set_point(const PointType& origin)
        {
          _origin = origin;
        }
      }; // class DistanceFunction
 
      template<int dim_, typename DataType_>
      class DistanceFunctionSD :
        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;
 
        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 _a, _b;
 
        public:
          explicit Evaluator(const DistanceFunctionSD& function) :
            _origin(function._origin),
            _a(function._a),
            _b(function._b)
          {
          }
 
          ValueType value(const PointType& point) const
          {
            return _a + _b*(point - _origin).norm_euclid();
          }
 
          GradientType gradient(const PointType& point) const
          {
            DataType norm ((point - _origin).norm_euclid());
            if(norm < Math::eps<DataType>())
              return GradientType::null();
 
            return (_b / norm) * (point - _origin);
          }
 
          HessianType hessian(const PointType& point) const
          {
            DataType norm ((point - _origin).norm_euclid());
            if(norm <= Math::eps<DataType>())
              return HessianType::null();
 
            HessianType hess(DataType(0));
 
            norm = DataType(1)/norm;
            DataType denom = _b*Math::sqr(norm)*norm;
 
            for(int i(0); i < dim_; ++i)
            {
              hess[i][i] = _b*norm;
              for(int j(0); j < dim_; ++j)
              {
                hess[i][j] -= ((point[i] - _origin[i]) * (point[j] - _origin[j]))*denom;
              }
            }
 
            return hess;
          }
        }; // class DistanceFunctionSD::Evaluator<...>
 
      private:
        PointType _origin;
        DataType_ _a;
        DataType_ _b;
 
      public:
        explicit DistanceFunctionSD(const PointType& origin, const DataType_ a_, const DataType_ b_) :
          _origin(origin),
          _a(a_),
          _b(b_)
        {
        }
 
        void set_point(const PointType& origin)
        {
          _origin = origin;
        }
      }; // class DistanceFunctionSD
 
      template<int component_, int dim_, typename DataType_>
      class PlaneDistanceFunctionSD :
        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;
 
        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_ _b;
 
        public:
          explicit Evaluator(const PlaneDistanceFunctionSD& function) :
            _origin(function._origin),
            _b(function._b)
          {
          }
 
          ValueType value(const PointType& point) const
          {
            return _b * Math::abs(point[component_] - _origin[component_]);
          }
 
          GradientType gradient(const PointType& point) const
          {
            DataType norm = this->value(point);
            if(norm < Math::eps<DataType>())
              return GradientType::null();
 
            GradientType grad;
            grad.format();
            grad[component_] = _b * Math::signum(point[component_] - _origin[component_]);
            return grad;
          }
 
          HessianType hessian(const PointType& DOXY(point)) const
          {
            return HessianType::null();
          }
        }; // class PlaneDistanceFunctionSD::Evaluator<...>
 
      private:
        PointType _origin;
        DataType_ _b;
 
      public:
        explicit PlaneDistanceFunctionSD(const PointType& origin, const DataType_ b_) :
          _origin(origin),
          _b(b_)
        {
        }
 
        void set_point(const PointType& origin)
        {
          _origin = origin;
        }
      }; // class PlaneDistanceFunctionSD
 
      template<int dim_, typename DataType_>
      class InverseDistanceFunction :
        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;
 
        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;
 
        public:
          explicit Evaluator(const InverseDistanceFunction& function) :
            _origin(function._origin)
          {
          }
 
          ValueType value(const PointType& point) const
          {
            return DataType(1) / (point - _origin).norm_euclid();
          }
 
          GradientType gradient(const PointType& point) const
          {
            const DataType norm_sqr = (point - _origin).norm_euclid_sqr();
 
            if(norm_sqr < Math::eps<DataType>())
              return GradientType::null();
 
            return (DataType(1) / Math::pow(norm_sqr, DataType(1.5))) * (_origin - point);
          }
 
          HessianType hessian(const PointType& point) const
          {
            const DataType norm_sqr = (point - _origin).norm_euclid_sqr();
            if(norm_sqr < Math::eps<DataType>())
              return HessianType::null();
 
            const DataType denom = DataType(1) / Math::pow(norm_sqr, DataType(2.5));
            HessianType hess (DataType(0));
            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 InverseDistanceFunction::Evaluator<...>
 
      public:
        PointType _origin;
 
      public:
        InverseDistanceFunction()
        {
          _origin.format();
        }
 
        explicit InverseDistanceFunction(const PointType& origin) :
          _origin(origin)
        {
        }
 
        void set_point(const PointType& origin)
        {
          _origin = origin;
        }
      }; // class InverseDistanceFunction
 
#ifdef FEAT_HAVE_CGAL
      template<typename DT_>
      class CGALDistanceFunction : public Analytic::Function
      {
      public:
        typedef DT_ DataType;
        typedef Analytic::Image::Scalar ImageType;
        static constexpr int domain_dim = 3;
        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;
 
      private:
        Geometry::CGALWrapper<DT_> _cgal_wrapper;
      public:
 
        explicit CGALDistanceFunction(std::istream& filestream_, Geometry::CGALFileMode filemode_) :
          _cgal_wrapper(filestream_, filemode_)
        {
        }
 
        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:
          typedef typename CGALDistanceFunction::DataType IDT_;
          const Geometry::CGALWrapper<IDT_>* _cgal_wrapper;
 
        public:
          explicit Evaluator(const CGALDistanceFunction& function) :
            _cgal_wrapper(&function._cgal_wrapper)
          {
          }
 
          ValueType value(const PointType& point)
          {
            return ((_cgal_wrapper->point_inside(IDT_(point[0]), IDT_(point[1]), IDT_(point[2])) ? ValueType(1) : ValueType(-1)) *
                    ValueType(Math::sqrt(_cgal_wrapper->squared_distance(IDT_(point[0]), IDT_(point[1]), IDT_(point[2])))));
          }
        }; // class CGALDistFunc::Evaluator<...>
      }; // class CGALDistFunc<...>
#endif // FEAT_HAVE_CGAL
    } // namespace Distance
  } // namespace Analytic
} // namespace FEAT