• Main Page
  • Namespaces
  • Data Structures
  • Files
  • File List
  • Globals

/export/development/ViennaGrid/release/ViennaGrid-1.0.0/viennagrid/point.hpp

Go to the documentation of this file.
00001 #ifndef VIENNAGRID_POINT_HPP
00002 #define VIENNAGRID_POINT_HPP
00003 
00004 /* =======================================================================
00005    Copyright (c) 2011, Institute for Microelectronics,
00006                        Institute for Analysis and Scientific Computing,
00007                        TU Wien.
00008 
00009                             -----------------
00010                      ViennaGrid - The Vienna Grid Library
00011                             -----------------
00012 
00013    Authors:      Karl Rupp                           rupp@iue.tuwien.ac.at
00014                  Josef Weinbub                    weinbub@iue.tuwien.ac.at
00015                
00016    (A list of additional contributors can be found in the PDF manual)
00017 
00018    License:      MIT (X11), see file LICENSE in the base directory
00019 ======================================================================= */
00020 
00021 
00022 #include <cmath>
00023 #include <assert.h>
00024 #include <stdexcept>
00025 #include <cstddef>
00026 #include <sstream>
00027 
00028 #include "viennagrid/forwards.h"
00029 #include "viennagrid/traits/point.hpp"
00030 
00035 namespace viennagrid
00036 {
00037   
00038   template <long d>
00039   struct dim_dispatcher;
00040   
00046   template <typename FromPointType, 
00047             typename ToPointType,
00048             typename FromCoordinateSystem = typename traits::coordinate_system<FromPointType>::type,
00049             typename ToCoordinateSystem = typename traits::coordinate_system<ToPointType>::type 
00050            >
00051   class coordinate_converter
00052   {
00053     public:
00055       ToPointType operator()(FromPointType const & p_in)
00056       {
00057         typedef typename FromPointType::ERROR_COORDINATE_SYSTEM_UNKNOWN   type;
00058       }
00059   };
00060   
00061   
00063   template <typename FromPointType, 
00064             typename ToPointType>
00065   class coordinate_converter<FromPointType, ToPointType, cartesian_cs<2>, polar_cs>
00066   {
00067     public:
00068       ToPointType operator()(FromPointType const & p_in)
00069       {
00070         ToPointType ret;
00071         ret[0] = sqrt(p_in[0] * p_in[0] + p_in[1] * p_in[1]);
00072         ret[1] = atan2(p_in[1], p_in[0]);
00073         return ret;
00074       }
00075   };
00076   
00078   template <typename FromPointType, 
00079             typename ToPointType>
00080   class coordinate_converter<FromPointType, ToPointType, cartesian_cs<3>, spherical_cs>
00081   {
00082     public:
00083       ToPointType operator()(FromPointType const & p_in)
00084       {
00085         ToPointType ret;
00086         ret[0] = sqrt(p_in[0] * p_in[0] + p_in[1] * p_in[1] + p_in[2] * p_in[2]);
00087         ret[1] = (ret[0] != 0) ? acos(p_in[2] / ret[0]) : 0;
00088         ret[2] = atan2(p_in[1], p_in[0]);
00089         return ret;
00090       }
00091   };
00092   
00094   template <typename FromPointType, 
00095             typename ToPointType>
00096   class coordinate_converter<FromPointType, ToPointType, cartesian_cs<3>, cylindrical_cs>
00097   {
00098     public:
00099       ToPointType operator()(FromPointType const & p_in)
00100       {
00101         ToPointType ret;
00102         ret[0] = sqrt(p_in[0] * p_in[0] + p_in[1] * p_in[1]);
00103         ret[1] = atan2(p_in[1], p_in[0]);
00104         ret[2] = p_in[2];
00105         return ret;
00106       }
00107   };
00108 
00109   namespace result_of
00110   {
00112     template <typename PointType>
00113     struct cartesian_point
00114     {
00115       typedef viennagrid::point_t<typename traits::value_type<PointType>::type,
00116                                   viennagrid::cartesian_cs<viennagrid::traits::static_size<PointType>::value>
00117                                  >                    type;
00118     };
00119 
00120   }
00121   
00123   template <typename PointType, typename CoordinateSystem>
00124   typename result_of::cartesian_point<PointType>::type
00125   to_cartesian_impl(PointType const & p, CoordinateSystem const &)
00126   {
00127     typedef typename result_of::cartesian_point<PointType>::type CartesianPointType;
00128     
00129     return coordinate_converter<PointType, CartesianPointType>()(p);
00130   }
00131 
00132   //public interface
00136   template <typename PointType>
00137   typename result_of::cartesian_point<PointType>::type
00138   to_cartesian(PointType const & p)
00139   {
00140     return to_cartesian_impl(p, typename traits::coordinate_system<PointType>::type());
00141   }
00142   
00143 
00145   template <typename FromPointType, 
00146             typename ToPointType>
00147   class coordinate_converter<FromPointType, ToPointType, polar_cs, cartesian_cs<2> >
00148   {
00149     public:
00150       ToPointType operator()(FromPointType const & p_in)
00151       {
00152         ToPointType ret;
00153         ret[0] = p_in[0] * cos(p_in[1]);
00154         ret[1] = p_in[0] * sin(p_in[1]);
00155         return ret;
00156       }
00157   };
00158   
00159 
00161   template <typename FromPointType, 
00162             typename ToPointType>
00163   class coordinate_converter<FromPointType, ToPointType, spherical_cs, cartesian_cs<3> >
00164   {
00165     public:
00166       ToPointType operator()(FromPointType const & p_in)
00167       {
00168         ToPointType ret;
00169         ret[0] = p_in[0] * sin(p_in[1]) * cos(p_in[2]);
00170         ret[1] = p_in[0] * sin(p_in[1]) * sin(p_in[2]);
00171         ret[2] = p_in[0] * cos(p_in[1]);
00172         return ret;
00173       }
00174   };
00175 
00177   template <typename FromPointType, 
00178             typename ToPointType>
00179   class coordinate_converter<FromPointType, ToPointType, spherical_cs, cylindrical_cs>
00180   {
00181     public:
00182       ToPointType operator()(FromPointType const & p_in)
00183       {
00184         ToPointType ret;
00185         ret[0] = p_in[0] * sin(p_in[1]);                 //rho
00186         ret[1] = p_in[2];                                //phi
00187         ret[2] = p_in[0] * cos(p_in[1]);                 //z
00188         return ret;
00189       }
00190   };
00191 
00192   
00194   template <typename FromPointType, 
00195             typename ToPointType>
00196   class coordinate_converter<FromPointType, ToPointType, cylindrical_cs, cartesian_cs<3> >
00197   {
00198     public:
00199       ToPointType operator()(FromPointType const & p_in)
00200       {
00201         ToPointType ret;
00202         ret[0] = p_in[0] * cos(p_in[1]);
00203         ret[1] = p_in[0] * sin(p_in[1]);
00204         ret[2] = p_in[2];
00205         return ret;
00206       }
00207   };
00208 
00210   template <typename FromPointType, 
00211             typename ToPointType>
00212   class coordinate_converter<FromPointType, ToPointType, cylindrical_cs, spherical_cs>
00213   {
00214     public:
00215       ToPointType operator()(FromPointType const & p_in)
00216       {
00217         ToPointType ret;
00218         ret[0] = sqrt(p_in[0] * p_in[0] + p_in[2] * p_in[2]);
00219         ret[1] = (ret[0] != 0) ? acos(p_in[2] / ret[0]) : 0;
00220         ret[2] = p_in[1];
00221         return ret;
00222       }
00223   };
00224   
00225   /********************* CoordinateSystem *****************/
00226   
00228   template <long d>
00229   struct cartesian_cs 
00230   {
00231     enum { dim = d };
00232     
00234     template <typename PointType>
00235     static PointType add(PointType const & p1, PointType const & p2)
00236     {
00237       assert(p1.size() == p2.size());
00238       
00239       PointType ret;
00240       typedef typename PointType::size_type      size_type;
00241       for (size_type i=0; i<ret.size(); ++i)
00242         ret[i] = p1[i] + p2[i];
00243       return ret;
00244     }
00245 
00247     template <typename PointType>
00248     static void inplace_add(PointType & p1, PointType const & p2)
00249     {
00250       assert(p1.size() == p2.size());
00251       
00252       typedef typename PointType::size_type      size_type;
00253       for (size_type i=0; i<p1.size(); ++i)
00254         p1[i] += p2[i];
00255     }
00256 
00258     template <typename PointType>
00259     static PointType subtract(PointType const & p1, PointType const & p2)
00260     {
00261       typedef typename PointType::size_type      size_type;
00262       assert(p1.size() == p2.size());
00263       
00264       PointType ret;
00265       for (size_type i=0; i<ret.size(); ++i)
00266         ret[i] = p1[i] - p2[i];
00267       return ret;
00268     }
00269 
00271     template <typename PointType>
00272     static void inplace_subtract(PointType & p1, PointType const & p2)
00273     {
00274       typedef typename PointType::size_type      size_type;
00275       assert(p1.size() == p2.size());
00276       
00277       for (size_type i=0; i<p1.size(); ++i)
00278         p1[i] -= p2[i];
00279     }
00280     
00282     template <typename PointType>
00283     static PointType & inplace_stretch(PointType & p1, typename traits::value_type<PointType>::type factor)
00284     {
00285       typedef typename PointType::size_type      size_type;
00286       for (size_type i=0; i<p1.size(); ++i)
00287         p1[i] *= factor;
00288       return p1;
00289     }
00290     
00291   };
00292   
00294   template <typename PointType, long d>
00295   PointType const &
00296   to_cartesian_impl(PointType const & p, cartesian_cs<d>)
00297   {
00298     return p;
00299   }
00300 
00302   template <typename PointType, long d>
00303   PointType &
00304   to_cartesian_impl(PointType & p, cartesian_cs<d>)
00305   {
00306     return p;
00307   }
00308 
00309   
00310   
00312   template <typename CSystem>
00313   struct cs_base
00314   {
00316     template <typename PointType>
00317     static PointType add(PointType const & p1, PointType const & p2)
00318     {
00319       static const int DIM = viennagrid::traits::static_size<PointType>::value;
00320       
00321       typedef point_t<typename PointType::value_type, 
00322                       cartesian_cs<DIM> 
00323                      >      CartesianPointType;
00324       
00325       assert(p1.size() == p2.size());
00326       
00327       coordinate_converter<PointType, CartesianPointType,
00328                           CSystem, cartesian_cs<DIM> >  cs_to_cartesian;
00329       coordinate_converter<CartesianPointType, PointType, 
00330                            cartesian_cs<DIM>, CSystem>  cartesian_to_cs;
00331       
00332       CartesianPointType cs1 = cs_to_cartesian(p1);
00333       CartesianPointType cs2 = cs_to_cartesian(p2);
00334       CartesianPointType cs_ret = cartesian_cs<DIM>::add(cs1, cs2);
00335       PointType ret = cartesian_to_cs(cs_ret);
00336       
00337       return ret;
00338     }
00339 
00341     template <typename PointType>
00342     static void inplace_add(PointType & p1, PointType const & p2)
00343     {
00344       static const int DIM = viennagrid::traits::static_size<PointType>::value;
00345       
00346       typedef point_t<typename PointType::value_type, 
00347                       cartesian_cs<DIM> 
00348                      >      CartesianPointType;
00349       
00350       assert(p1.size() == p2.size());
00351       
00352       coordinate_converter<PointType, CartesianPointType, CSystem, cartesian_cs<DIM> >  cs_to_cartesian;
00353       coordinate_converter<CartesianPointType, PointType, cartesian_cs<DIM>, CSystem>  cartesian_to_cs;
00354       
00355       CartesianPointType cs1 = cs_to_cartesian(p1);
00356       CartesianPointType cs2 = cs_to_cartesian(p2);
00357       cartesian_cs<DIM>::inplace_add(cs1, cs2);
00358       p1 = cartesian_to_cs(cs1);
00359     }
00360 
00362     template <typename PointType>
00363     static PointType subtract(PointType const & p1, PointType const & p2)
00364     {
00365       static const int DIM = viennagrid::traits::static_size<PointType>::value;
00366       
00367       typedef point_t<typename PointType::value_type, 
00368                       cartesian_cs<DIM> 
00369                      >      CartesianPointType;
00370       
00371       assert(p1.size() == p2.size());
00372       
00373       coordinate_converter<PointType, CartesianPointType, CSystem, cartesian_cs<DIM> >  cs_to_cartesian;
00374       coordinate_converter<CartesianPointType, PointType, cartesian_cs<DIM>, CSystem>  cartesian_to_cs;
00375       
00376       CartesianPointType cs1 = cs_to_cartesian(p1);
00377       CartesianPointType cs2 = cs_to_cartesian(p2);
00378       CartesianPointType cs_ret = cartesian_cs<DIM>::subtract(cs1, cs2);
00379       PointType ret = cartesian_to_cs(cs_ret);
00380       
00381       return ret;
00382     }
00383 
00385     template <typename PointType>
00386     static void inplace_subtract(PointType & p1, PointType const & p2)
00387     {
00388       static const int DIM = viennagrid::traits::static_size<PointType>::value;
00389       
00390       typedef point_t<typename PointType::value_type, 
00391                       cartesian_cs<DIM> 
00392                      >      CartesianPointType;
00393       
00394       assert(p1.size() == p2.size());
00395       
00396       coordinate_converter<PointType, CartesianPointType, CSystem, cartesian_cs<DIM> >  cs_to_cartesian;
00397       coordinate_converter<CartesianPointType, PointType, cartesian_cs<DIM>, CSystem>  cartesian_to_cs;
00398       
00399       CartesianPointType cs1 = cs_to_cartesian(p1);
00400       CartesianPointType cs2 = cs_to_cartesian(p2);
00401       cartesian_cs<DIM>::inplace_subtract(cs1, cs2);
00402       p1 = cartesian_to_cs(cs1);
00403     }
00404   };
00405   
00407   struct polar_cs : public cs_base<polar_cs> 
00408   {
00409     enum { dim = 2 };
00410     
00411     template <typename PointType>
00412     static PointType & inplace_stretch(PointType & p1, typename traits::value_type<PointType>::type factor)
00413     {
00414       p1[0] *= factor;
00415       return p1;
00416     }
00417   };
00418   
00420   struct spherical_cs : public cs_base<spherical_cs> 
00421   {
00422     enum { dim = 3 };
00423     
00424     template <typename PointType>
00425     static PointType & inplace_stretch(PointType & p1, typename traits::value_type<PointType>::type factor)
00426     {
00427       p1[0] *= factor;
00428       return p1;
00429     }
00430   };
00431 
00433   struct cylindrical_cs : public cs_base<cylindrical_cs> 
00434   {
00435     enum { dim = 3 };
00436     
00438     template <typename PointType>
00439     static PointType & inplace_stretch(PointType & p1, typename traits::value_type<PointType>::type factor)
00440     {
00441       static const int DIM = viennagrid::traits::static_size<PointType>::value;
00442       
00443       typedef point_t<typename PointType::value_type, 
00444                       cartesian_cs<DIM> 
00445                      >      CartesianPointType;
00446       
00447       coordinate_converter<PointType, CartesianPointType, cylindrical_cs, cartesian_cs<DIM> >  cs_to_cartesian;
00448       coordinate_converter<CartesianPointType, PointType, cartesian_cs<DIM>, cylindrical_cs>  cartesian_to_cs;
00449       
00450       CartesianPointType cs1 = cs_to_cartesian(p1);
00451       cartesian_cs<DIM>::inplace_stretch(cs1, factor);
00452       p1 = cartesian_to_cs(cs1);
00453       return p1;
00454     }
00455   };
00456   
00457   
00458   
00459   
00460   /***************************** Point Type ************************/
00461   
00463   class point_index_out_of_bounds_exception : public std::exception
00464   {
00465     public:
00466       point_index_out_of_bounds_exception(std::size_t i) : i_(i) {};
00467       
00468       virtual const char* what() const throw()
00469       {
00470         std::stringstream ss;
00471         ss << "Point index " << i_ << " out of bounds!";
00472         return ss.str().c_str();
00473       }
00474       
00475     private:
00476       std::size_t i_;
00477   };
00478   
00479 
00481   template <typename CoordType, long d>
00482   struct point_filler
00483   {
00484     static void apply(CoordType * coords, CoordType x, CoordType y, CoordType z)
00485     {
00486       coords[0] = x;
00487       coords[1] = y;
00488       coords[2] = z;
00489       for (long i=3; i<d; ++i)
00490         coords[i] = 0;
00491     }
00492   };
00493   
00495   template <typename CoordType>
00496   struct point_filler<CoordType, 1>
00497   {
00498     static void apply(CoordType * coords, CoordType x, CoordType, CoordType)
00499     {
00500       coords[0] = x;
00501     }
00502   };
00503   
00505   template <typename CoordType>
00506   struct point_filler<CoordType, 2>
00507   {
00508     static void apply(CoordType * coords, CoordType x, CoordType y, CoordType)
00509     {
00510       coords[0] = x;
00511       coords[1] = y;
00512     }
00513   };
00514   
00516   template <typename CoordType>
00517   struct point_filler<CoordType, 3>
00518   {
00519     static void apply(CoordType * coords, CoordType x, CoordType y, CoordType z)
00520     {
00521       coords[0] = x;
00522       coords[1] = y;
00523       coords[2] = z;
00524     }
00525   };
00526   
00527   
00528   
00536   template <typename CoordType, typename CoordinateSystem>
00537   class point_t
00538   {
00539     public:
00541       typedef CoordType       value_type;
00543       typedef dim_type        size_type;
00544       
00546       enum { dim = CoordinateSystem::dim };
00547       
00549       point_t() 
00550       {
00551         point_filler<CoordType, dim>::apply(coords, 0, 0, 0);  //make sure that there is no bogus in the coords-array
00552       }
00553       
00555       point_t(CoordType x, CoordType y = 0, CoordType z = 0)
00556       {
00557         point_filler<CoordType, dim>::apply(coords, x, y, z);
00558       }
00559       
00561       template <typename CoordType2, typename CoordinateSystem2>
00562       point_t(point_t<CoordType2, CoordinateSystem2> const & p2)
00563       {
00564         //std::cout << "Copy CTOR!" << std::endl;
00565         *this = coordinate_converter<point_t<CoordType2, CoordinateSystem2>, point_t>()(p2);
00566       }
00567 
00568       //explicit copy CTOR
00570       point_t(point_t const & other)
00571       {
00572         for (size_type i=0; i<size_type(dim); ++i)
00573           coords[i] = other.coords[i];
00574       }
00575 
00577       point_t & operator=(point_t const & p2)
00578       {
00579         for (size_type i=0; i<size_type(dim); ++i)
00580           coords[i] = p2.coords[i];
00581         return *this;
00582       }
00583       
00584       
00586       template <typename CoordType2, typename CoordinateSystem2>
00587       point_t & operator=(point_t<CoordType2, CoordinateSystem2> const & p2)
00588       {
00589         //std::cout << "Assignment operator!" << std::endl;
00590         *this = coordinate_converter<point_t<CoordType2, CoordinateSystem2>, point_t>()(p2);
00591         return *this;
00592       }
00593       
00595       CoordType & operator[](size_type index) 
00596       {
00597         assert(index < dim);
00598         return coords[index]; 
00599       }
00600       
00602       CoordType const & operator[](size_type index) const
00603       {
00604         assert(index < dim);
00605         return coords[index]; 
00606       }
00607       
00609       CoordType at(size_type index)
00610       {
00611         if (index < 0 || index >= dim)
00612           throw point_index_out_of_bounds_exception(index);
00613         
00614         return coords[index];
00615       }
00616 
00618       CoordType const & at(size_type index) const
00619       {
00620         if (index < 0 || index >= dim)
00621           throw point_index_out_of_bounds_exception(index);
00622         
00623         return coords[index];
00624       }
00625       
00627       size_type size() const { return dim; }
00628       
00629       //
00630       // operators:
00631       //
00632       
00633       //with point:
00635       point_t operator+(point_t const & other) const
00636       {
00637         return CoordinateSystem::add(*this, other);
00638       }
00639 
00641       point_t & operator+=(point_t const & other)
00642       {
00643         CoordinateSystem::inplace_add(*this, other);
00644         return *this;
00645       }
00646 
00648       point_t operator-(point_t const & other) const
00649       {
00650         return CoordinateSystem::subtract(*this, other);
00651       }
00652       
00654       point_t & operator-=(point_t const & other)
00655       {
00656         CoordinateSystem::inplace_subtract(*this, other);
00657         return *this;
00658       }
00659       
00660       
00661       //with CoordType
00663       point_t & operator*=(CoordType factor)
00664       {
00665         CoordinateSystem::inplace_stretch(*this, factor);
00666         //for (size_type i=0; i<d; ++i)
00667         //  coords[i] *= factor;
00668         return *this;
00669       }
00670       
00672       point_t & operator/=(CoordType factor)
00673       {
00674         CoordinateSystem::inplace_stretch(*this, 1.0 / factor);
00675         //for (size_type i=0; i<d; ++i)
00676         //  coords[i] /= factor;
00677         return *this;
00678       }
00679 
00681       point_t operator*(CoordType factor) const
00682       {
00683         point_t ret(*this);
00684         return CoordinateSystem::inplace_stretch(ret, factor);
00685         //for (size_type i=0; i<d; ++i)
00686         //  ret[i] = coords[i] * factor;
00687         //return ret;
00688       }
00689 
00691       point_t operator/(CoordType factor) const
00692       {
00693         point_t ret(*this);
00694         return CoordinateSystem::inplace_stretch(ret, 1.0 / factor);
00695         //point_t ret;
00696         //for (size_type i=0; i<d; ++i)
00697         //  ret[i] = coords[i] / factor;
00698         //return ret;
00699       }
00700 
00701     private:
00702       CoordType coords[dim];
00703   };
00704 
00706   template <typename CoordType, typename CoordinateSystem>
00707   point_t<CoordType, CoordinateSystem>
00708   operator*(double val, point_t<CoordType, CoordinateSystem> const & p)
00709   {
00710     return p * val;
00711   }
00712 
00713 /*  template <typename CoordType, long d, typename CoordinateSystem>
00714   point_t<CoordType, d, CoordinateSystem>
00715   operator*(CoordType val, point_t<CoordType, d, CoordinateSystem> const & p)
00716   {
00717     return p * val;
00718   }*/
00719 
00721   template <typename CoordType, typename CoordinateSystem>
00722   std::ostream& operator << (std::ostream & os, point_t<CoordType, CoordinateSystem> const & p)
00723   {
00724     typedef typename point_t<CoordType, CoordinateSystem>::size_type      size_type;
00725     for (size_type i=0; i<static_cast<size_type>(CoordinateSystem::dim); ++i)
00726       os << p[i] << " ";
00727     return os;
00728   }
00729   
00730 
00732   struct point_less
00733   {
00734     template <typename PointType>
00735     bool operator()(PointType const & p1, PointType const & p2) const
00736     {
00737       for (std::size_t i=0; i<p1.size(); ++i)
00738       {
00739         if (p1[i] < p2[i])
00740           return true;
00741         else if (p1[i] > p2[i])
00742           return false;
00743       }
00744       return false;
00745     }
00746   };
00747 
00748 }
00749 #endif

Generated on Wed Sep 14 2011 19:21:30 for ViennaGrid by  doxygen 1.7.1