viennagrid/point.hpp

Go to the documentation of this file.

#ifndef VIENNAGRID_POINT_HPP
#define VIENNAGRID_POINT_HPP

/* =======================================================================
   Copyright (c) 2011-2014, Institute for Microelectronics,
                            Institute for Analysis and Scientific Computing,
                            TU Wien.

                            -----------------
                     ViennaGrid - The Vienna Grid Library
                            -----------------

   License:      MIT (X11), see file LICENSE in the base directory
======================================================================= */


#include <cmath>
#include <assert.h>
#include <stdexcept>
#include <cstddef>
#include <sstream>

#include "viennagrid/forwards.hpp"

#include "viennagrid/storage/static_array.hpp"

namespace viennagrid
{
  namespace result_of
  {
    template <typename PointType>
    struct dimension;

    template <typename CoordType, typename CoordinateSystem>
    struct dimension< spatial_point<CoordType, CoordinateSystem> >
    {
      static const int value = CoordinateSystem::dim;
    };


    template <typename PointType>
    struct coordinate_system
    {
      //by default, we don't know anything about the point type, so let's complain at compile time
      typedef typename PointType::ERROR_UNKNOWN_COORDINATE_SYSTEM_FOR_POINT_TYPE    type;
    };

    template <typename CoordType, typename CoordinateSystem>
    struct coordinate_system< spatial_point<CoordType, CoordinateSystem> >
    {
      typedef CoordinateSystem    type;
    };


    template <typename PointType>
    struct static_size;

    template <typename CoordType, typename CoordinateSystem>
    struct static_size< spatial_point<CoordType, CoordinateSystem> >
    {
      static const int value = CoordinateSystem::dim;
    };

  }

  template <typename PointType>
  std::size_t dynamic_size(PointType const & p)
  {
    return p.size();
  }

  template <int d>
  struct dim_dispatcher;

  template <typename FromPointType,
            typename ToPointType,
            typename FromCoordinateSystem = typename result_of::coordinate_system<FromPointType>::type,
            typename ToCoordinateSystem = typename result_of::coordinate_system<ToPointType>::type
           >
  class coordinate_converter
  {
  public:
    ToPointType operator()(FromPointType const &)
    {
      typename FromPointType::ERROR_COORDINATE_SYSTEM_UNKNOWN   error_object;
      (void)error_object;
    }
  };


  template <typename FromPointType,
            typename ToPointType>
  class coordinate_converter<FromPointType, ToPointType, cartesian_cs<2>, polar_cs>
  {
  public:
    ToPointType operator()(FromPointType const & p_in)
    {
      ToPointType ret;
      ret[0] = sqrt(p_in[0] * p_in[0] + p_in[1] * p_in[1]);
      ret[1] = atan2(p_in[1], p_in[0]);
      return ret;
    }
  };

  template <typename FromPointType,
            typename ToPointType>
  class coordinate_converter<FromPointType, ToPointType, cartesian_cs<3>, spherical_cs>
  {
  public:
    ToPointType operator()(FromPointType const & p_in)
    {
      ToPointType ret;
      ret[0] = sqrt(p_in[0] * p_in[0] + p_in[1] * p_in[1] + p_in[2] * p_in[2]);
      ret[1] = (std::fabs(ret[0]) > 0) ? acos(p_in[2] / ret[0]) : 0;
      ret[2] = atan2(p_in[1], p_in[0]);
      return ret;
    }
  };

  template <typename FromPointType,
            typename ToPointType>
  class coordinate_converter<FromPointType, ToPointType, cartesian_cs<3>, cylindrical_cs>
  {
  public:
    ToPointType operator()(FromPointType const & p_in)
    {
      ToPointType ret;
      ret[0] = sqrt(p_in[0] * p_in[0] + p_in[1] * p_in[1]);
      ret[1] = atan2(p_in[1], p_in[0]);
      ret[2] = p_in[2];
      return ret;
    }
  };

  namespace result_of
  {
    template <typename PointType>
    struct cartesian_point
    {
      typedef viennagrid::spatial_point<typename result_of::coord<PointType>::type,
                                  viennagrid::cartesian_cs<viennagrid::result_of::static_size<PointType>::value>
                                 >                    type;
    };

  }

  template <typename PointType, typename CoordinateSystem>
  typename result_of::cartesian_point<PointType>::type
  to_cartesian_impl(PointType const & p, CoordinateSystem const &)
  {
    typedef typename result_of::cartesian_point<PointType>::type CartesianPointType;

    return coordinate_converter<PointType, CartesianPointType>()(p);
  }


  template <typename FromPointType,
            typename ToPointType>
  class coordinate_converter<FromPointType, ToPointType, polar_cs, cartesian_cs<2> >
  {
  public:
    ToPointType operator()(FromPointType const & p_in)
    {
      ToPointType ret;
      ret[0] = p_in[0] * cos(p_in[1]);
      ret[1] = p_in[0] * sin(p_in[1]);
      return ret;
    }
  };


  template <typename FromPointType,
            typename ToPointType>
  class coordinate_converter<FromPointType, ToPointType, spherical_cs, cartesian_cs<3> >
  {
  public:
    ToPointType operator()(FromPointType const & p_in)
    {
      ToPointType ret;
      ret[0] = p_in[0] * sin(p_in[1]) * cos(p_in[2]);
      ret[1] = p_in[0] * sin(p_in[1]) * sin(p_in[2]);
      ret[2] = p_in[0] * cos(p_in[1]);
      return ret;
    }
  };

  template <typename FromPointType,
            typename ToPointType>
  class coordinate_converter<FromPointType, ToPointType, spherical_cs, cylindrical_cs>
  {
  public:
    ToPointType operator()(FromPointType const & p_in)
    {
      ToPointType ret;
      ret[0] = p_in[0] * sin(p_in[1]);                 //rho
      ret[1] = p_in[2];                                //phi
      ret[2] = p_in[0] * cos(p_in[1]);                 //z
      return ret;
    }
  };


  template <typename FromPointType,
            typename ToPointType>
  class coordinate_converter<FromPointType, ToPointType, cylindrical_cs, cartesian_cs<3> >
  {
  public:
    ToPointType operator()(FromPointType const & p_in)
    {
      ToPointType ret;
      ret[0] = p_in[0] * cos(p_in[1]);
      ret[1] = p_in[0] * sin(p_in[1]);
      ret[2] = p_in[2];
      return ret;
    }
  };

  template <typename FromPointType,
            typename ToPointType>
  class coordinate_converter<FromPointType, ToPointType, cylindrical_cs, spherical_cs>
  {
  public:
    ToPointType operator()(FromPointType const & p_in)
    {
      ToPointType ret;
      ret[0] = sqrt(p_in[0] * p_in[0] + p_in[2] * p_in[2]);
      ret[1] = (std::fabs(ret[0]) > 0) ? acos(p_in[2] / ret[0]) : 0;
      ret[2] = p_in[1];
      return ret;
    }
  };

  /********************* CoordinateSystem *****************/

  template <int d>
  struct cartesian_cs
  {
    static const int dim = d;

    template <typename PointType>
    static PointType add(PointType const & p1, PointType const & p2)
    {
      assert(p1.size() == p2.size());

      PointType ret;
      typedef typename PointType::size_type      size_type;
      for (size_type i=0; i<ret.size(); ++i)
        ret[i] = p1[i] + p2[i];
      return ret;
    }

    template <typename PointType>
    static void inplace_add(PointType & p1, PointType const & p2)
    {
      assert(p1.size() == p2.size());

      typedef typename PointType::size_type      size_type;
      for (size_type i=0; i<p1.size(); ++i)
        p1[i] += p2[i];
    }

    template <typename PointType>
    static PointType subtract(PointType const & p1, PointType const & p2)
    {
      typedef typename PointType::size_type      size_type;
      assert(p1.size() == p2.size());

      PointType ret;
      for (size_type i=0; i<ret.size(); ++i)
        ret[i] = p1[i] - p2[i];
      return ret;
    }

    template <typename PointType>
    static void inplace_subtract(PointType & p1, PointType const & p2)
    {
      typedef typename PointType::size_type      size_type;
      assert(p1.size() == p2.size());

      for (size_type i=0; i<p1.size(); ++i)
        p1[i] -= p2[i];
    }

    template <typename PointType>
    static PointType & inplace_stretch(PointType & p1, typename result_of::coord<PointType>::type factor)
    {
      typedef typename PointType::size_type      size_type;
      for (size_type i=0; i<p1.size(); ++i)
        p1[i] *= factor;
      return p1;
    }

  };

  template <typename PointType, int d>
  PointType const &
  to_cartesian_impl(PointType const & p, cartesian_cs<d>)
  {
    return p;
  }

  template <typename PointType, int d>
  PointType &
  to_cartesian_impl(PointType & p, cartesian_cs<d>)
  {
    return p;
  }


  //public interface
  template <typename PointType>
  typename result_of::cartesian_point<PointType>::type
  to_cartesian(PointType const & p)
  {
    return to_cartesian_impl(p, typename result_of::coordinate_system<PointType>::type());
  }



  template <typename CSystem>
  struct cs_base
  {
    template <typename PointType>
    static PointType add(PointType const & p1, PointType const & p2)
    {
      static const int DIM = viennagrid::result_of::static_size<PointType>::value;

      typedef spatial_point<typename PointType::value_type,
                      cartesian_cs<DIM>
                     >      CartesianPointType;

      assert(p1.size() == p2.size());

      coordinate_converter<PointType, CartesianPointType,
                          CSystem, cartesian_cs<DIM> >  cs_to_cartesian;
      coordinate_converter<CartesianPointType, PointType,
                           cartesian_cs<DIM>, CSystem>  cartesian_to_cs;

      CartesianPointType cs1 = cs_to_cartesian(p1);
      CartesianPointType cs2 = cs_to_cartesian(p2);
      CartesianPointType cs_ret = cartesian_cs<DIM>::add(cs1, cs2);
      PointType ret = cartesian_to_cs(cs_ret);

      return ret;
    }

    template <typename PointType>
    static void inplace_add(PointType & p1, PointType const & p2)
    {
      static const int DIM = viennagrid::result_of::static_size<PointType>::value;

      typedef spatial_point<typename PointType::value_type,
                      cartesian_cs<DIM>
                     >      CartesianPointType;

      assert(p1.size() == p2.size());

      coordinate_converter<PointType, CartesianPointType, CSystem, cartesian_cs<DIM> >  cs_to_cartesian;
      coordinate_converter<CartesianPointType, PointType, cartesian_cs<DIM>, CSystem>  cartesian_to_cs;

      CartesianPointType cs1 = cs_to_cartesian(p1);
      CartesianPointType cs2 = cs_to_cartesian(p2);
      cartesian_cs<DIM>::inplace_add(cs1, cs2);
      p1 = cartesian_to_cs(cs1);
    }

    template <typename PointType>
    static PointType subtract(PointType const & p1, PointType const & p2)
    {
      static const int DIM = viennagrid::result_of::static_size<PointType>::value;

      typedef spatial_point<typename PointType::value_type,
                      cartesian_cs<DIM>
                     >      CartesianPointType;

      assert(p1.size() == p2.size());

      coordinate_converter<PointType, CartesianPointType, CSystem, cartesian_cs<DIM> >  cs_to_cartesian;
      coordinate_converter<CartesianPointType, PointType, cartesian_cs<DIM>, CSystem>  cartesian_to_cs;

      CartesianPointType cs1 = cs_to_cartesian(p1);
      CartesianPointType cs2 = cs_to_cartesian(p2);
      CartesianPointType cs_ret = cartesian_cs<DIM>::subtract(cs1, cs2);
      PointType ret = cartesian_to_cs(cs_ret);

      return ret;
    }

    template <typename PointType>
    static void inplace_subtract(PointType & p1, PointType const & p2)
    {
      static const int DIM = viennagrid::result_of::static_size<PointType>::value;

      typedef spatial_point<typename PointType::value_type,
                      cartesian_cs<DIM>
                     >      CartesianPointType;

      assert(p1.size() == p2.size());

      coordinate_converter<PointType, CartesianPointType, CSystem, cartesian_cs<DIM> >  cs_to_cartesian;
      coordinate_converter<CartesianPointType, PointType, cartesian_cs<DIM>, CSystem>  cartesian_to_cs;

      CartesianPointType cs1 = cs_to_cartesian(p1);
      CartesianPointType cs2 = cs_to_cartesian(p2);
      cartesian_cs<DIM>::inplace_subtract(cs1, cs2);
      p1 = cartesian_to_cs(cs1);
    }
  };

  struct polar_cs : public cs_base<polar_cs>
  {
    static const int dim = 2;

    template <typename PointType>
    static PointType & inplace_stretch(PointType & p1, typename result_of::coord<PointType>::type factor)
    {
      p1[0] *= factor;
      return p1;
    }
  };

  struct spherical_cs : public cs_base<spherical_cs>
  {
    static const int dim = 3;

    template <typename PointType>
    static PointType & inplace_stretch(PointType & p1, typename result_of::coord<PointType>::type factor)
    {
      p1[0] *= factor;
      return p1;
    }
  };

  struct cylindrical_cs : public cs_base<cylindrical_cs>
  {
    static const int dim = 3;

    template <typename PointType>
    static PointType & inplace_stretch(PointType & p1, typename result_of::coord<PointType>::type factor)
    {
      static const int DIM = viennagrid::result_of::static_size<PointType>::value;

      typedef spatial_point<typename PointType::value_type,
                      cartesian_cs<DIM>
                     >      CartesianPointType;

      coordinate_converter<PointType, CartesianPointType, cylindrical_cs, cartesian_cs<DIM> >  cs_to_cartesian;
      coordinate_converter<CartesianPointType, PointType, cartesian_cs<DIM>, cylindrical_cs>  cartesian_to_cs;

      CartesianPointType cs1 = cs_to_cartesian(p1);
      cartesian_cs<DIM>::inplace_stretch(cs1, factor);
      p1 = cartesian_to_cs(cs1);
      return p1;
    }
  };




  /***************************** Point Type ************************/

  class point_index_out_of_bounds_exception : public std::exception
  {
  public:
    point_index_out_of_bounds_exception(std::size_t i) : i_(i) {}

    virtual const char* what() const throw()
    {
      std::stringstream ss;
      ss << "Point index " << i_ << " out of bounds!";
      return ss.str().c_str();
    }

  private:
    std::size_t i_;
  };


  template <typename CoordType, int d>
  struct point_filler
  {
    static void apply(CoordType * coords, CoordType x, CoordType y, CoordType z)
    {
      coords[0] = x;
      coords[1] = y;
      coords[2] = z;
      for (int i=3; i<d; ++i)
        coords[i] = 0;
    }
  };

  template <typename CoordType>
  struct point_filler<CoordType, 1>
  {
    static void apply(CoordType * coords, CoordType x, CoordType, CoordType)
    {
      coords[0] = x;
    }
  };

  template <typename CoordType>
  struct point_filler<CoordType, 2>
  {
    static void apply(CoordType * coords, CoordType x, CoordType y, CoordType)
    {
      coords[0] = x;
      coords[1] = y;
    }
  };

  template <typename CoordType>
  struct point_filler<CoordType, 3>
  {
    static void apply(CoordType * coords, CoordType x, CoordType y, CoordType z)
    {
      coords[0] = x;
      coords[1] = y;
      coords[2] = z;
    }
  };



  template <typename CoordType, typename CoordinateSystem>
  class spatial_point : public static_array<CoordType, CoordinateSystem::dim>
  {
  public:
    typedef CoordType       value_type;
    typedef dim_type        size_type;


    static const int dim = CoordinateSystem::dim;

    spatial_point()
    {
      point_filler<CoordType, dim>::apply( &(*this)[0], 0, 0, 0);  //make sure that there is no bogus in the coords-array
    }

    spatial_point(CoordType x, CoordType y = 0, CoordType z = 0)
    {
      point_filler<CoordType, dim>::apply( &(*this)[0], x, y, z);
    }

    template <typename CoordType2, typename CoordinateSystem2>
    spatial_point(spatial_point<CoordType2, CoordinateSystem2> const & p2)
    {
      *this = coordinate_converter<spatial_point<CoordType2, CoordinateSystem2>, spatial_point>()(p2);
    }


    template <typename CoordType2, typename CoordinateSystem2>
    spatial_point & operator=(spatial_point<CoordType2, CoordinateSystem2> const & p2)
    {
      *this = coordinate_converter<spatial_point<CoordType2, CoordinateSystem2>, spatial_point>()(p2);
      return *this;
    }

    //
    // operators:
    //

    spatial_point operator-() const
    {
      return spatial_point() - *this;
    }

    //with point:
    spatial_point operator+(spatial_point const & other) const
    {
      return CoordinateSystem::add(*this, other);
    }

    spatial_point & operator+=(spatial_point const & other)
    {
      CoordinateSystem::inplace_add(*this, other);
      return *this;
    }

    spatial_point operator-(spatial_point const & other) const
    {
      return CoordinateSystem::subtract(*this, other);
    }

    spatial_point & operator-=(spatial_point const & other)
    {
      CoordinateSystem::inplace_subtract(*this, other);
      return *this;
    }


    //with CoordType
    spatial_point & operator*=(CoordType factor)
    {
      CoordinateSystem::inplace_stretch(*this, factor);
      return *this;
    }

    spatial_point & operator/=(CoordType factor)
    {
      CoordinateSystem::inplace_stretch(*this, 1.0 / factor);
      return *this;
    }

    spatial_point operator*(CoordType factor) const
    {
      spatial_point ret(*this);
      return CoordinateSystem::inplace_stretch(ret, factor);
    }

    spatial_point operator/(CoordType factor) const
    {
      spatial_point ret(*this);
      return CoordinateSystem::inplace_stretch(ret, 1.0 / factor);
    }

  };


  template<typename CoordType, typename CoordinateSystem>
  spatial_point<CoordType, CoordinateSystem> min(const spatial_point<CoordType, CoordinateSystem> & p1, const spatial_point<CoordType, CoordinateSystem> & p2)
  {
    spatial_point<CoordType, CoordinateSystem> tmp;
    for (std::size_t i = 0; i < spatial_point<CoordType, CoordinateSystem>::dim; ++i)
      tmp[i] = std::min(p1[i], p2[i]);
    return tmp;
  }

  template<typename CoordType, typename CoordinateSystem>
  spatial_point<CoordType, CoordinateSystem> max(const spatial_point<CoordType, CoordinateSystem> & p1, const spatial_point<CoordType, CoordinateSystem> & p2)
  {
    spatial_point<CoordType, CoordinateSystem> tmp;
    for (std::size_t i = 0; i < spatial_point<CoordType, CoordinateSystem>::dim; ++i)
      tmp[i] = std::max(p1[i], p2[i]);
    return tmp;
  }


  template <typename CoordType, typename CoordinateSystem>
  spatial_point<CoordType, CoordinateSystem>
  operator*(double val, spatial_point<CoordType, CoordinateSystem> const & p)
  {
    return p * val;
  }

  template <typename CoordType, typename CoordinateSystem>
  std::ostream& operator << (std::ostream & os, spatial_point<CoordType, CoordinateSystem> const & p)
  {
    typedef typename spatial_point<CoordType, CoordinateSystem>::size_type      size_type;
    os << "(";
    for (size_type i=0; i<static_cast<size_type>(CoordinateSystem::dim); ++i)
      os << p[i] << (i == static_cast<size_type>(CoordinateSystem::dim)-1 ? "" :" ");
    os << ")";
    return os;
  }


  struct point_less
  {
    template <typename PointType>
    bool operator()(PointType const & p1, PointType const & p2) const
    {
      for (std::size_t i=0; i<p1.size(); ++i)
      {
        if (p1[i] < p2[i])
          return true;
        else if (p1[i] > p2[i])
          return false;
      }
      return false;
    }
  };

  namespace result_of
  {
    template<typename CoordType, typename CoordinateSystem>
    struct point< spatial_point<CoordType, CoordinateSystem> >
    {
      typedef spatial_point<CoordType, CoordinateSystem> type;
    };

    template<typename CoordType, typename CoordinateSystem>
    struct coord< spatial_point<CoordType, CoordinateSystem> >
    {
      typedef CoordType type;
    };

    template<typename CoordType, typename CoordinateSystem>
    struct geometric_dimension< spatial_point<CoordType, CoordinateSystem> >
    {
      static const int value = spatial_point<CoordType, CoordinateSystem>::dim;
    };
  }

}
#endif