Rheolef  7.1
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 //#include "rheolef/distributed.h"
82 
83 namespace rheolef {
84 
85 // [verbatim_point_basic]
86 template <class T>
87 class point_basic {
88 public:
89 
90 // typedefs:
91 
92  typedef size_t size_type;
93  typedef T element_type;
94  typedef T scalar_type;
95  typedef T float_type;
96 
97 // allocators:
98 
99  explicit point_basic();
100  explicit point_basic (const T& x0, const T& x1 = 0, const T& x2 = 0);
101 
102  template <class T1>
104 
105  template <class T1>
107 
108  point_basic (const std::initializer_list<T>& il);
109 
110 // accessors:
111 
112  T& operator[](int i_coord) { return _x[i_coord%3]; }
113  T& operator()(int i_coord) { return _x[i_coord%3]; }
114  const T& operator[](int i_coord) const { return _x[i_coord%3]; }
115  const T& operator()(int i_coord) const { return _x[i_coord%3]; }
116 
117 // algebra:
118 
119  bool operator== (const point_basic<T>& v) const;
120  bool operator!= (const point_basic<T>& v) const;
128 
129  template <class U>
130  typename
131  std::enable_if<
134  >::type
135  operator* (const U& a) const;
136  point_basic<T> operator/ (const T& a) const;
138 
139 // i/o:
140 
141  std::istream& get (std::istream& s, int d = 3);
142  std::ostream& put (std::ostream& s, int d = 3) const;
143 // [verbatim_point_basic]
144 
145 // interface for CGAL library inter-operability:
146 
147  const T& x() const { return _x[0]; }
148  const T& y() const { return _x[1]; }
149  const T& z() const { return _x[2]; }
150  T& x(){ return _x[0]; }
151  T& y(){ return _x[1]; }
152  T& z(){ return _x[2]; }
153 
154 // data:
155 // protected:
156  T _x[3];
157 // internal:
158  static T _my_abs(const T& x) { return (x > T(0)) ? x : -x; }
159 // [verbatim_point_basic_cont]
160 };
161 // [verbatim_point_basic_cont]
162 
163 // [verbatim_point]
165 // [verbatim_point]
166 
167 // [verbatim_point_basic_cont2]
168 template<class T>
169 std::istream& operator >> (std::istream& s, point_basic<T>& p);
170 
171 template<class T>
172 std::ostream& operator << (std::ostream& s, const point_basic<T>& p);
173 
174 template <class T, class U>
175 typename
176 std::enable_if<
177  details::is_rheolef_arithmetic<U>::value
179 >::type
180 operator* (const U& a, const point_basic<T>& u);
181 
182 template<class T>
184 vect (const point_basic<T>& v, const point_basic<T>& w);
185 
186 // metrics:
187 template<class T>
188 T dot (const point_basic<T>& x, const point_basic<T>& y);
189 
190 template<class T>
191 T norm2 (const point_basic<T>& x);
192 
193 template<class T>
194 T norm (const point_basic<T>& x);
195 
196 template<class T>
197 T dist2 (const point_basic<T>& x, const point_basic<T>& y);
198 
199 template<class T>
200 T dist (const point_basic<T>& x, const point_basic<T>& y);
201 
202 template<class T>
203 T dist_infty (const point_basic<T>& x, const point_basic<T>& y);
204 
205 template <class T>
206 T vect2d (const point_basic<T>& v, const point_basic<T>& w);
207 
208 template <class T>
209 T mixt (const point_basic<T>& u, const point_basic<T>& v, const point_basic<T>& w);
210 
211 // robust(exact) floating point predicates: return the sign of the value as (0, > 0, < 0)
212 // formally: orient2d(a,b,x) = vect2d(a-x,b-x)
213 template <class T>
214 int sign_orient2d (
215  const point_basic<T>& a,
216  const point_basic<T>& b,
217  const point_basic<T>& c);
218 
219 template <class T>
220 int sign_orient3d (
221  const point_basic<T>& a,
222  const point_basic<T>& b,
223  const point_basic<T>& c,
224  const point_basic<T>& d);
225 
226 // compute also the value:
227 template <class T>
228 T orient2d(
229  const point_basic<T>& a,
230  const point_basic<T>& b,
231  const point_basic<T>& c);
232 
233 // formally: orient3d(a,b,c,x) = mixt3d(a-x,b-x,c-x)
234 template <class T>
235 T orient3d(
236  const point_basic<T>& a,
237  const point_basic<T>& b,
238  const point_basic<T>& c,
239  const point_basic<T>& d);
240 
241 template <class T>
242 std::string ptos (const point_basic<T>& x, int d = 3);
243 
244 // ccomparators: lexicographic order
245 template<class T, size_t d>
247 // [verbatim_point_basic_cont2]
248 // -------------------------------------------------------------------------------------
249 // inline'd
250 // -------------------------------------------------------------------------------------
251 template <class T, class U>
252 inline
253 typename
254 std::enable_if<
255  details::is_rheolef_arithmetic<U>::value
257 >::type
258 operator* (const U& a, const point_basic<T>& u)
259 {
260  return point_basic<T> (a*u[0], a*u[1], a*u[2]);
261 }
262 template<class T>
263 inline
265 vect (const point_basic<T>& v, const point_basic<T>& w)
266 {
267  return point_basic<T> (
268  v[1]*w[2]-v[2]*w[1],
269  v[2]*w[0]-v[0]*w[2],
270  v[0]*w[1]-v[1]*w[0]);
271 }
272 // metrics:
273 template<class T>
274 inline
275 T dot (const point_basic<T>& x, const point_basic<T>& y)
276 {
277  return x[0]*y[0]+x[1]*y[1]+x[2]*y[2];
278 }
279 template<class T>
280 inline
282 {
283  return dot(x,x);
284 }
285 template<class T>
286 inline
288 {
289  return sqrt(norm2(x));
290 }
291 template<class T>
292 inline
293 T dist2 (const point_basic<T>& x, const point_basic<T>& y)
294 {
295  return norm2(x-y);
296 }
297 template<class T>
298 inline
299 T dist (const point_basic<T>& x, const point_basic<T>& y)
300 {
301  return norm(x-y);
302 }
303 template<class T>
304 inline
306 {
307  return max(point_basic<T>::_my_abs(x[0]-y[0]),
308  max(point_basic<T>::_my_abs(x[1]-y[1]),
309  point_basic<T>::_my_abs(x[2]-y[2])));
310 }
311 // ccomparators: lexicographic order
312 template<class T, size_t d>
313 inline
314 bool
316 {
317  for (typename point_basic<T>::size_type i = 0; i < d; i++) {
318  if (a[i] < b[i]) return true;
319  if (a[i] > b[i]) return false;
320  }
321  return false; // equality
322 }
324 template<class T> struct scalar_traits { typedef T type; };
325 template<class T> struct scalar_traits<point_basic<T> > { typedef T type; };
326 template<class T> struct float_traits<point_basic<T> > { typedef typename float_traits<T>::type type; };
327 
328 template<class T>
330  _x[0] = T();
331  _x[1] = T();
332  _x[2] = T();
333 }
334 template<class T>
336  const T& x0,
337  const T& x1,
338  const T& x2)
339 {
340  _x[0] = x0;
341  _x[1] = x1;
342  _x[2] = x2;
343 }
344 template<class T>
345 template<class T1>
347 {
348  _x[0] = p._x[0];
349  _x[1] = p._x[1];
350  _x[2] = p._x[2];
351 }
352 template<class T>
353 template <class T1>
356 {
357  _x[0] = p._x[0];
358  _x[1] = p._x[1];
359  _x[2] = p._x[2];
360  return *this;
361 }
362 template<class T>
363 point_basic<T>::point_basic (const std::initializer_list<T>& il) : _x() {
364  typedef typename std::initializer_list<T>::const_iterator const_iterator;
365  check_macro (il.size() <= 3, "unexpected initializer list size=" << il.size() << " > 3");
366  size_type i = 0;
367  for (const_iterator iter = il.begin(); iter != il.end(); ++iter, ++i) {
368  _x[i] = *iter;
369  }
370  for (i = il.size(); i < 3; ++i) {
371  _x[i] = T();
372  }
373 }
374 // input/output
375 template<class T>
376 std::istream&
377 point_basic<T>::get (std::istream& s, int d)
378 {
379  switch (d) {
380  case 0 : _x[0] = _x[1] = _x[2] = T(0); return s;
381  case 1 : _x[1] = _x[2] = T(0); return s >> _x[0];
382  case 2 : _x[2] = T(0); return s >> _x[0] >> _x[1];
383  default: return s >> _x[0] >> _x[1] >> _x[2];
384  }
385 }
386 template<class T>
387 inline
388 std::ostream&
389 point_basic<T>::put (std::ostream& s, int d) const
390 {
391  switch (d) {
392  case 0 : return s;
393  case 1 : return s << _x[0];
394  case 2 : return s << _x[0] << " " << _x[1];
395  default: return s << _x[0] << " " << _x[1] << " " << _x[2];
396  }
397 }
398 template<class T>
399 inline
400 std::istream&
401 operator >> (std::istream& s, point_basic<T>& p)
402 {
403  return s >> p[0] >> p[1] >> p[2];
404 }
405 template<class T>
406 inline
407 std::ostream&
408 operator << (std::ostream& s, const point_basic<T>& p)
409 {
410  return s << p[0] << " " << p[1] << " " << p[2];
411 }
412 template<class T>
413 std::string
414 ptos (const point_basic<T>& x, int d)
415 {
416  std::ostringstream ostrstr;
417  x.put (ostrstr, d);
418  return ostrstr.str();
419 }
420 // ----------------------------------------------------------
421 // point extra: inlined
422 // ----------------------------------------------------------
423 #define def_point_function2(f,F) \
424 template<class T> \
425 inline \
426 point_basic<T> \
427 f (const point_basic<T>& x) \
428 { \
429  point_basic<T> y; \
430  for (size_t i = 0; i < 3; i++) \
431  y[i] = F(x[i]); \
432  return y; \
433 }
434 
435 #define def_point_function(f) def_point_function2(f,f)
436 
438 def_point_function(sqrt)
440 def_point_function(log10)
442 #undef def_point_function2
443 #undef def_point_function
444 
445 template <class T>
446 bool
448 {
449  return _x[0] == v[0] && _x[1] == v[1] && _x[2] == v[2];
450 }
451 template <class T>
452 bool
454 {
455  return !operator==(v);
456 }
457 template <class T>
460 {
461  _x[0] += v[0];
462  _x[1] += v[1];
463  _x[2] += v[2];
464  return *this;
465 }
466 template <class T>
469 {
470  _x[0] -= v[0];
471  _x[1] -= v[1];
472  _x[2] -= v[2];
473  return *this;
474 }
475 template <class T>
478 {
479  _x[0] *= a;
480  _x[1] *= a;
481  _x[2] *= a;
482  return *this;
483 }
484 template <class T>
487 {
488  _x[0] /= a;
489  _x[1] /= a;
490  _x[2] /= a;
491  return *this;
492 }
493 template <class T>
496 {
497  return point_basic<T> (_x[0]+v[0],
498  _x[1]+v[1],
499  _x[2]+v[2]);
500 }
501 template <class T>
504 {
505  return point_basic<T> (-_x[0],
506  -_x[1],
507  -_x[2]);
508 }
509 template <class T>
512 {
513  return point_basic<T> (_x[0]-v[0],
514  _x[1]-v[1],
515  _x[2]-v[2]);
516 }
517 template<class T1, class T2>
518 inline
520 operator/ (const T2& a, const point_basic<T1>& x)
521 {
522  point_basic<T1> y;
523  for (size_t i = 0; i < 3; i++)
524  if (x[i] != 0) y[i] = a/x[i];
525  return y;
526 }
527 template <class T>
528 template <class U>
529 typename
530 std::enable_if<
531  details::is_rheolef_arithmetic<U>::value
533 >::type
534 point_basic<T>::operator* (const U& a) const
535 {
536  return point_basic<T> (_x[0]*a,
537  _x[1]*a,
538  _x[2]*a);
539 }
540 template <class T>
543 {
544  return operator* (T(1)/T(a));
545 }
546 template <class T>
549 {
550  return point_basic<T> (_x[0]/v[0],
551  _x[1]/v[1],
552  _x[2]/v[2]);
553 }
554 // vect2d() and mixt() are deduced from:
555 template <class T>
556 inline
557 T
559 {
560  return orient2d (point_basic<T>(), v, w);
561 }
562 template <class T>
563 inline
564 T
565 mixt (const point_basic<T>& u, const point_basic<T>& v, const point_basic<T>& w)
566 {
567  return orient3d (point_basic<T>(), u, v, w);
568 }
569 
570 }// namespace rheolef
571 // -------------------------------------------------------------
572 // serialization
573 // -------------------------------------------------------------
574 #ifdef _RHEOLEF_HAVE_MPI
575 #include <boost/serialization/serialization.hpp>
576 
577 namespace boost {
578  namespace serialization {
579  template<class Archive, class T>
580  void serialize (Archive& ar, class rheolef::point_basic<T>& x, const unsigned int version) {
581  ar & x[0];
582  ar & x[1];
583  ar & x[2];
584  }
585  } // namespace serialization
586 } // namespace boost
587 
588 // Some serializable types, like geo_element, have a fixed amount of data stored at fixed field positions.
589 // When this is the case, boost::mpi can optimize their serialization and transmission to avoid extraneous copy operations.
590 // To enable this optimization, we specialize the type trait is_mpi_datatype, e.g.:
591 namespace boost {
592  namespace mpi {
593  // TODO: when point_basic<T> is not a simple type, such as T=bigfloat or T=gmp, etc
594  template <>
595  struct is_mpi_datatype<rheolef::point_basic<double> > : mpl::true_ { };
596  } // namespace mpi
597 } // namespace boost
598 #endif // _RHEOLEF_HAVE_MPI
599 
600 #endif /* _RHEO_BASIC_POINT_H */
point_basic< T > operator-() const
Definition: point.h:503
const T & y() const
Definition: point.h:148
std::istream & get(std::istream &s, int d=3)
Definition: point.h:377
bool operator==(const point_basic< T > &v) const
const T & operator[](int i_coord) const
Definition: point.h:114
size_t size_type
Definition: point.h:92
point_basic< T > & operator+=(const point_basic< T > &v)
point_basic< T > operator+(const point_basic< T > &v) const
Definition: point.h:495
T & operator[](int i_coord)
Definition: point.h:112
std::ostream & put(std::ostream &s, int d=3) const
Definition: point.h:389
T & operator()(int i_coord)
Definition: point.h:113
bool operator!=(const point_basic< T > &v) const
point_basic< T > operator/(const T &a) const
Definition: point.h:542
const T & z() const
Definition: point.h:149
point_basic< T > & operator=(const point_basic< T1 > &p)
Definition: point.h:355
static T _my_abs(const T &x)
Definition: point.h:158
point_basic< T > & operator*=(const T &a)
Definition: point.h:477
point_basic< T > & operator/=(const T &a)
Definition: point.h:486
const T & x() const
Definition: point.h:147
point_basic(const T &x0, const T &x1=0, const T &x2=0)
Definition: point.h:335
point_basic(const std::initializer_list< T > &il)
Definition: point.h:363
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:115
point_basic< T > & operator-=(const point_basic< T > &v)
Definition: point.h:468
point_basic< Float > point
Definition: point.h:164
rheolef::std type
point_basic< T >
Definition: piola_fem.h:135
Expr1::float_type T
Definition: field_expr.h:261
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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:299
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:558
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:305
T mixt(const point_basic< T > &u, const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:565
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition: point.h:293
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:315
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:265
std::string ptos(const point_basic< T > &x, int d=3)
Definition: point.h:414
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:437
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:326
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:324
Definition: leveque.h:25