Rheolef  7.2
an efficient C++ finite element environment
point.h
Go to the documentation of this file.
1 # ifndef _RHEO_BASIC_POINT_H
2 # define _RHEO_BASIC_POINT_H
3 //
4 // This file is part of Rheolef.
5 //
6 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
7 //
8 // Rheolef is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // Rheolef is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with Rheolef; if not, write to the Free Software
20 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // =========================================================================
23 // author: Pierre.Saramito@imag.fr
24 
25 namespace rheolef {
77 } // namespace rheolef
78 
79 #include "rheolef/compiler_mpi.h"
80 #include <sstream>
81 
82 namespace rheolef {
83 
84 // [verbatim_point_basic]
85 template <class T>
86 class point_basic {
87 public:
88 
89 // typedefs:
90 
91  typedef size_t size_type;
92  typedef T element_type;
93  typedef T scalar_type;
94  typedef T float_type;
95 
96 // allocators:
97 
98  explicit point_basic();
99  explicit point_basic (const T& x0, const T& x1 = 0, const T& x2 = 0);
100 
101  template <class T1>
103 
104  template <class T1>
106 
107  point_basic (const std::initializer_list<T>& il);
108 
109 // accessors:
110 
111  T& operator[](int i_coord) { return _x[i_coord%3]; }
112  T& operator()(int i_coord) { return _x[i_coord%3]; }
113  const T& operator[](int i_coord) const { return _x[i_coord%3]; }
114  const T& operator()(int i_coord) const { return _x[i_coord%3]; }
115 
116 // algebra:
117 
118  bool operator== (const point_basic<T>& v) const;
119  bool operator!= (const point_basic<T>& v) const;
127 
128  template <class U>
129  typename
130  std::enable_if<
133  >::type
134  operator* (const U& a) const;
135  point_basic<T> operator/ (const T& a) const;
137 
138 // i/o:
139 
140  std::istream& get (std::istream& s, int d = 3);
141  std::ostream& put (std::ostream& s, int d = 3) const;
142 // [verbatim_point_basic]
143 
144 // interface for CGAL library inter-operability:
145 
146  const T& x() const { return _x[0]; }
147  const T& y() const { return _x[1]; }
148  const T& z() const { return _x[2]; }
149  T& x(){ return _x[0]; }
150  T& y(){ return _x[1]; }
151  T& z(){ return _x[2]; }
152 
153 // data:
154 // protected:
155  T _x[3];
156 // internal:
157  static T _my_abs(const T& x) { return (x > T(0)) ? x : -x; }
158 // [verbatim_point_basic_cont]
159 };
160 // [verbatim_point_basic_cont]
161 
162 // [verbatim_point]
164 // [verbatim_point]
165 
166 // [verbatim_point_basic_cont2]
167 template<class T>
168 std::istream& operator >> (std::istream& s, point_basic<T>& p);
169 
170 template<class T>
171 std::ostream& operator << (std::ostream& s, const point_basic<T>& p);
172 
173 template <class T, class U>
174 typename
175 std::enable_if<
176  details::is_rheolef_arithmetic<U>::value
178 >::type
179 operator* (const U& a, const point_basic<T>& u);
180 
181 template<class T>
183 vect (const point_basic<T>& v, const point_basic<T>& w);
184 
185 // metrics:
186 template<class T>
187 T dot (const point_basic<T>& x, const point_basic<T>& y);
188 
189 template<class T>
190 T norm2 (const point_basic<T>& x);
191 
192 template<class T>
193 T norm (const point_basic<T>& x);
194 
195 template<class T>
196 T dist2 (const point_basic<T>& x, const point_basic<T>& y);
197 
198 template<class T>
199 T dist (const point_basic<T>& x, const point_basic<T>& y);
200 
201 template<class T>
202 T dist_infty (const point_basic<T>& x, const point_basic<T>& y);
203 
204 template <class T>
205 T vect2d (const point_basic<T>& v, const point_basic<T>& w);
206 
207 template <class T>
208 T mixt (const point_basic<T>& u, const point_basic<T>& v, const point_basic<T>& w);
209 
210 // robust(exact) floating point predicates: return the sign of the value as (0, > 0, < 0)
211 // formally: orient2d(a,b,x) = vect2d(a-x,b-x)
212 template <class T>
213 int sign_orient2d (
214  const point_basic<T>& a,
215  const point_basic<T>& b,
216  const point_basic<T>& c);
217 
218 template <class T>
219 int sign_orient3d (
220  const point_basic<T>& a,
221  const point_basic<T>& b,
222  const point_basic<T>& c,
223  const point_basic<T>& d);
224 
225 // compute also the value:
226 template <class T>
227 T orient2d(
228  const point_basic<T>& a,
229  const point_basic<T>& b,
230  const point_basic<T>& c);
231 
232 // formally: orient3d(a,b,c,x) = mixt3d(a-x,b-x,c-x)
233 template <class T>
234 T orient3d(
235  const point_basic<T>& a,
236  const point_basic<T>& b,
237  const point_basic<T>& c,
238  const point_basic<T>& d);
239 
240 template <class T>
241 std::string ptos (const point_basic<T>& x, int d = 3);
242 
243 // ccomparators: lexicographic order
244 template<class T, size_t d>
246 // [verbatim_point_basic_cont2]
247 // -------------------------------------------------------------------------------------
248 // inline'd
249 // -------------------------------------------------------------------------------------
250 template <class T, class U>
251 inline
252 typename
253 std::enable_if<
254  details::is_rheolef_arithmetic<U>::value
256 >::type
257 operator* (const U& a, const point_basic<T>& u)
258 {
259  return point_basic<T> (a*u[0], a*u[1], a*u[2]);
260 }
261 template<class T>
262 inline
264 vect (const point_basic<T>& v, const point_basic<T>& w)
265 {
266  return point_basic<T> (
267  v[1]*w[2]-v[2]*w[1],
268  v[2]*w[0]-v[0]*w[2],
269  v[0]*w[1]-v[1]*w[0]);
270 }
271 // metrics:
272 template<class T>
273 inline
274 T dot (const point_basic<T>& x, const point_basic<T>& y)
275 {
276  return x[0]*y[0]+x[1]*y[1]+x[2]*y[2];
277 }
278 template<class T>
279 inline
281 {
282  return dot(x,x);
283 }
284 template<class T>
285 inline
287 {
288  return sqrt(norm2(x));
289 }
290 template<class T>
291 inline
292 T dist2 (const point_basic<T>& x, const point_basic<T>& y)
293 {
294  return norm2(x-y);
295 }
296 template<class T>
297 inline
298 T dist (const point_basic<T>& x, const point_basic<T>& y)
299 {
300  return norm(x-y);
301 }
302 template<class T>
303 inline
305 {
306  return max(point_basic<T>::_my_abs(x[0]-y[0]),
307  max(point_basic<T>::_my_abs(x[1]-y[1]),
308  point_basic<T>::_my_abs(x[2]-y[2])));
309 }
310 // ccomparators: lexicographic order
311 template<class T, size_t d>
312 inline
313 bool
315 {
316  for (typename point_basic<T>::size_type i = 0; i < d; i++) {
317  if (a[i] < b[i]) return true;
318  if (a[i] > b[i]) return false;
319  }
320  return false; // equality
321 }
323 template<class T> struct scalar_traits { typedef T type; };
324 template<class T> struct scalar_traits<point_basic<T> > { typedef T type; };
325 template<class T> struct float_traits<point_basic<T> > { typedef typename float_traits<T>::type type; };
326 
327 template<class T>
329  _x[0] = T();
330  _x[1] = T();
331  _x[2] = T();
332 }
333 template<class T>
335  const T& x0,
336  const T& x1,
337  const T& x2)
338 {
339  _x[0] = x0;
340  _x[1] = x1;
341  _x[2] = x2;
342 }
343 template<class T>
344 template<class T1>
346 {
347  _x[0] = p._x[0];
348  _x[1] = p._x[1];
349  _x[2] = p._x[2];
350 }
351 template<class T>
352 template <class T1>
355 {
356  _x[0] = p._x[0];
357  _x[1] = p._x[1];
358  _x[2] = p._x[2];
359  return *this;
360 }
361 template<class T>
362 point_basic<T>::point_basic (const std::initializer_list<T>& il) : _x() {
363  typedef typename std::initializer_list<T>::const_iterator const_iterator;
364  check_macro (il.size() <= 3, "unexpected initializer list size=" << il.size() << " > 3");
365  size_type i = 0;
366  for (const_iterator iter = il.begin(); iter != il.end(); ++iter, ++i) {
367  _x[i] = *iter;
368  }
369  for (i = il.size(); i < 3; ++i) {
370  _x[i] = T();
371  }
372 }
373 // input/output
374 template<class T>
375 std::istream&
376 point_basic<T>::get (std::istream& s, int d)
377 {
378  switch (d) {
379  case 0 : _x[0] = _x[1] = _x[2] = T(0); return s;
380  case 1 : _x[1] = _x[2] = T(0); return s >> _x[0];
381  case 2 : _x[2] = T(0); return s >> _x[0] >> _x[1];
382  default: return s >> _x[0] >> _x[1] >> _x[2];
383  }
384 }
385 template<class T>
386 inline
387 std::ostream&
388 point_basic<T>::put (std::ostream& s, int d) const
389 {
390  switch (d) {
391  case 0 : return s;
392  case 1 : return s << _x[0];
393  case 2 : return s << _x[0] << " " << _x[1];
394  default: return s << _x[0] << " " << _x[1] << " " << _x[2];
395  }
396 }
397 template<class T>
398 inline
399 std::istream&
400 operator >> (std::istream& s, point_basic<T>& p)
401 {
402  return s >> p[0] >> p[1] >> p[2];
403 }
404 template<class T>
405 inline
406 std::ostream&
407 operator << (std::ostream& s, const point_basic<T>& p)
408 {
409  return s << p[0] << " " << p[1] << " " << p[2];
410 }
411 template<class T>
412 std::string
413 ptos (const point_basic<T>& x, int d)
414 {
415  std::ostringstream ostrstr;
416  x.put (ostrstr, d);
417  return ostrstr.str();
418 }
419 // ----------------------------------------------------------
420 // point extra: inlined
421 // ----------------------------------------------------------
422 #define def_point_function2(f,F) \
423 template<class T> \
424 inline \
425 point_basic<T> \
426 f (const point_basic<T>& x) \
427 { \
428  point_basic<T> y; \
429  for (size_t i = 0; i < 3; i++) \
430  y[i] = F(x[i]); \
431  return y; \
432 }
433 
434 #define def_point_function(f) def_point_function2(f,f)
435 
437 def_point_function(sqrt)
439 def_point_function(log10)
441 #undef def_point_function2
442 #undef def_point_function
443 
444 template <class T>
445 bool
447 {
448  return _x[0] == v[0] && _x[1] == v[1] && _x[2] == v[2];
449 }
450 template <class T>
451 bool
453 {
454  return !operator==(v);
455 }
456 template <class T>
459 {
460  _x[0] += v[0];
461  _x[1] += v[1];
462  _x[2] += v[2];
463  return *this;
464 }
465 template <class T>
468 {
469  _x[0] -= v[0];
470  _x[1] -= v[1];
471  _x[2] -= v[2];
472  return *this;
473 }
474 template <class T>
477 {
478  _x[0] *= a;
479  _x[1] *= a;
480  _x[2] *= a;
481  return *this;
482 }
483 template <class T>
486 {
487  _x[0] /= a;
488  _x[1] /= a;
489  _x[2] /= a;
490  return *this;
491 }
492 template <class T>
495 {
496  return point_basic<T> (_x[0]+v[0],
497  _x[1]+v[1],
498  _x[2]+v[2]);
499 }
500 template <class T>
503 {
504  return point_basic<T> (-_x[0],
505  -_x[1],
506  -_x[2]);
507 }
508 template <class T>
511 {
512  return point_basic<T> (_x[0]-v[0],
513  _x[1]-v[1],
514  _x[2]-v[2]);
515 }
516 template<class T1, class T2>
517 inline
519 operator/ (const T2& a, const point_basic<T1>& x)
520 {
521  point_basic<T1> y;
522  for (size_t i = 0; i < 3; i++)
523  if (x[i] != 0) y[i] = a/x[i];
524  return y;
525 }
526 template <class T>
527 template <class U>
528 typename
529 std::enable_if<
530  details::is_rheolef_arithmetic<U>::value
532 >::type
533 point_basic<T>::operator* (const U& a) const
534 {
535  return point_basic<T> (_x[0]*a,
536  _x[1]*a,
537  _x[2]*a);
538 }
539 template <class T>
542 {
543  return operator* (T(1)/T(a));
544 }
545 template <class T>
548 {
549  return point_basic<T> (_x[0]/v[0],
550  _x[1]/v[1],
551  _x[2]/v[2]);
552 }
553 // vect2d() and mixt() are deduced from:
554 template <class T>
555 inline
556 T
558 {
559  return orient2d (point_basic<T>(), v, w);
560 }
561 template <class T>
562 inline
563 T
564 mixt (const point_basic<T>& u, const point_basic<T>& v, const point_basic<T>& w)
565 {
566  return orient3d (point_basic<T>(), u, v, w);
567 }
568 
569 }// namespace rheolef
570 // -------------------------------------------------------------
571 // serialization
572 // -------------------------------------------------------------
573 #ifdef _RHEOLEF_HAVE_MPI
574 #include <boost/serialization/serialization.hpp>
575 
576 namespace boost {
577  namespace serialization {
578  template<class Archive, class T>
579  void serialize (Archive& ar, class rheolef::point_basic<T>& x, const unsigned int version) {
580  ar & x[0];
581  ar & x[1];
582  ar & x[2];
583  }
584  } // namespace serialization
585 } // namespace boost
586 
587 // Some serializable types, like geo_element, have a fixed amount of data stored at fixed field positions.
588 // When this is the case, boost::mpi can optimize their serialization and transmission to avoid extraneous copy operations.
589 // To enable this optimization, we specialize the type trait is_mpi_datatype, e.g.:
590 namespace boost {
591  namespace mpi {
592  // TODO: when point_basic<T> is not a simple type, such as T=bigfloat or T=gmp, etc
593  template <>
594  struct is_mpi_datatype<rheolef::point_basic<double> > : mpl::true_ { };
595  } // namespace mpi
596 } // namespace boost
597 #endif // _RHEOLEF_HAVE_MPI
598 
599 #endif /* _RHEO_BASIC_POINT_H */
point_basic< T > operator-() const
Definition: point.h:502
const T & y() const
Definition: point.h:147
std::istream & get(std::istream &s, int d=3)
Definition: point.h:376
bool operator==(const point_basic< T > &v) const
const T & operator[](int i_coord) const
Definition: point.h:113
size_t size_type
Definition: point.h:91
point_basic< T > & operator+=(const point_basic< T > &v)
point_basic< T > operator+(const point_basic< T > &v) const
Definition: point.h:494
T & operator[](int i_coord)
Definition: point.h:111
std::ostream & put(std::ostream &s, int d=3) const
Definition: point.h:388
T & operator()(int i_coord)
Definition: point.h:112
bool operator!=(const point_basic< T > &v) const
point_basic< T > operator/(const T &a) const
Definition: point.h:541
const T & z() const
Definition: point.h:148
point_basic< T > & operator=(const point_basic< T1 > &p)
Definition: point.h:354
static T _my_abs(const T &x)
Definition: point.h:157
point_basic< T > & operator*=(const T &a)
Definition: point.h:476
point_basic< T > & operator/=(const T &a)
Definition: point.h:485
const T & x() const
Definition: point.h:146
point_basic(const T &x0, const T &x1=0, const T &x2=0)
Definition: point.h:334
point_basic(const std::initializer_list< T > &il)
Definition: point.h:362
std::enable_if< details::is_rheolef_arithmetic< U >::value,point_basic< T > >::type operator*(const U &a) const
const T & operator()(int i_coord) const
Definition: point.h:114
point_basic< T > & operator-=(const point_basic< T > &v)
Definition: point.h:467
point_basic< Float > point
Definition: point.h:163
rheolef::std type
point_basic< T >
Definition: piola_fem.h:135
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void serialize(Archive &ar, class rheolef::point_basic< T > &x, const unsigned int version)
Definition: point.h:579
This file is part of Rheolef.
bool operator==(const heap_allocator< T1 > &lhs, const heap_allocator< T1 > &rhs)
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition: point.h:298
std::istream & operator>>(std::istream &is, const catchmark &m)
Definition: catchmark.h:88
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
T vect2d(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:557
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
T dist_infty(const point_basic< T > &x, const point_basic< T > &y)
Definition: point.h:304
T mixt(const point_basic< T > &u, const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:564
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition: point.h:292
dia< T, M > operator/(const T &lambda, const dia< T, M > &d)
Definition: dia.h:145
T orient2d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c)
bool lexicographically_less(const point_basic< T > &a, const point_basic< T > &b)
Definition: point.h:314
T orient3d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c, const point_basic< T > &d)
int sign_orient3d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c, const point_basic< T > &d)
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:264
std::string ptos(const point_basic< T > &x, int d=3)
Definition: point.h:413
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
int sign_orient2d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c)
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
def_point_function(sqr) def_point_function(sqrt) def_point_function(log) def_point_function(log10) def_point_function(exp) template< class T > bool point_basic< T >
Definition: point.h:436
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition: tensor-exp.cc:92
Definition: sphere.icc:25
float_traits< T >::type type
Definition: point.h:325
helper for std::complex<T>: get basic T type
Definition: Float.h:93
helper for point_basic<T> & tensor_basic<T>: get basic T type
Definition: point.h:323
Definition: leveque.h:25