00001 #ifndef VIENNAGRID_POINT_HPP
00002 #define VIENNAGRID_POINT_HPP
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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
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]);
00186 ret[1] = p_in[2];
00187 ret[2] = p_in[0] * cos(p_in[1]);
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
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
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);
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
00565 *this = coordinate_converter<point_t<CoordType2, CoordinateSystem2>, point_t>()(p2);
00566 }
00567
00568
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
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
00631
00632
00633
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
00663 point_t & operator*=(CoordType factor)
00664 {
00665 CoordinateSystem::inplace_stretch(*this, factor);
00666
00667
00668 return *this;
00669 }
00670
00672 point_t & operator/=(CoordType factor)
00673 {
00674 CoordinateSystem::inplace_stretch(*this, 1.0 / factor);
00675
00676
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
00686
00687
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
00696
00697
00698
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
00714
00715
00716
00717
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