Rheolef  7.1
an efficient C++ finite element environment
point_predicate.cc
Go to the documentation of this file.
1 // robust floating point predicates:
21 // implementation using exact CGAL predicates, when available
22 // together witha custom cgal kernel that uses rheolef::point_basic<T>
23 //
24 // author: Pierre.Saramito@imag.fr
25 //
26 // date: 12 march 2012
27 //
28 // all the work for exact predicates is performed by :
29 // "CGAL/internal/Static_filters/Orientation_3.h"
30 //
31 #include "rheolef/point.h"
32 
33 #ifdef _RHEOLEF_HAVE_CGAL
34 #include "rheolef/cgal_traits.h"
35 #endif // _RHEOLEF_HAVE_CGAL
36 
37 namespace rheolef {
38 
39 template<class T>
40 static
41 T
42 inexact_orient2d (const point_basic<T>& x, const point_basic<T>& a,
43  const point_basic<T>& b)
44 {
45  T ax0 = a[0] - x[0];
46  T bx0 = b[0] - x[0];
47  T ax1 = a[1] - x[1];
48  T bx1 = b[1] - x[1];
49  return ax0*bx1 - ax1*bx0;
50 }
51 template <class T>
52 static
53 T
54 inexact_orient3d (const point_basic<T>& x, const point_basic<T>& a,
55  const point_basic<T>& b, const point_basic<T>& c)
56 {
57  T ax0 = a[0] - x[0];
58  T bx0 = b[0] - x[0];
59  T cx0 = c[0] - x[0];
60  T ax1 = a[1] - x[1];
61  T bx1 = b[1] - x[1];
62  T cx1 = c[1] - x[1];
63  T ax2 = a[2] - x[2];
64  T bx2 = b[2] - x[2];
65  T cx2 = c[2] - x[2];
66  return ax0 * (bx1 * cx2 - bx2 * cx1)
67  + bx0 * (cx1 * ax2 - cx2 * ax1)
68  + cx0 * (ax1 * bx2 - ax2 * bx1);
69 }
70 template <class T>
71 int
73  const point_basic<T>& a,
74  const point_basic<T>& b,
75  const point_basic<T>& c)
76 {
77 #ifdef _RHEOLEF_HAVE_CGAL
78  typedef typename geo_cgal_traits<T,2>::Kernel Kernel;
79  typename Kernel::Orientation_2 orientation;
80  CGAL::Orientation sgn = orientation(a, b, c);
81  return (sgn == CGAL::NEGATIVE) ? -1 : ((sgn == CGAL::ZERO) ? 0 : 1);
82 #else // _RHEOLEF_HAVE_CGAL
83  T sgn = inexact_orient2d(a, b, c);
84  return (sgn < 0) ? -1 : ((sgn == 0) ? 0 : 1);
85 #endif // _RHEOLEF_HAVE_CGAL
86 }
87 template <class T>
88 int
90  const point_basic<T>& a,
91  const point_basic<T>& b,
92  const point_basic<T>& c,
93  const point_basic<T>& d)
94 {
95 #ifdef _RHEOLEF_HAVE_CGAL
96  typedef typename geo_cgal_traits<T,3>::Kernel Kernel;
97  typename Kernel::Orientation_3 orientation;
98  CGAL::Orientation sgn = orientation(a, b, c, d);
99  return (sgn == CGAL::NEGATIVE) ? -1 : ((sgn == CGAL::ZERO) ? 0 : 1);
100 #else // _RHEOLEF_HAVE_CGAL
101  T sgn = inexact_orient3d(a, b, c, d);
102  return (sgn < 0) ? -1 : ((sgn == 0) ? 0 : 1);
103 #endif // _RHEOLEF_HAVE_CGAL
104 }
105 template<class T>
106 T
108  const point_basic<T>& c)
109 {
110  T value = inexact_orient2d(a, b, c);
111 #ifndef _RHEOLEF_HAVE_CGAL
112  return value;
113 #else // _RHEOLEF_HAVE_CGAL
114  int sgn = sign_orient2d(a, b, c);
115  value = fabs(value);
116  if (sgn == 0) return 0;
117  if (value != 0) return sgn*value;
118  // sgn != 0 but value == 0
120 #endif // _RHEOLEF_HAVE_CGAL
121 }
122 template <class T>
123 T
125  const point_basic<T>& c, const point_basic<T>& d)
126 {
127  T value = inexact_orient3d(a, b, c, d);
128 #ifndef _RHEOLEF_HAVE_CGAL
129  return value;
130 #else // _RHEOLEF_HAVE_CGAL
131  int sgn = sign_orient3d(a, b, c, d);
132  value = fabs(value);
133  if (sgn == 0) return 0;
134  if (value != 0) return sgn*value;
135  // sgn != 0 but value == 0
137 #endif // _RHEOLEF_HAVE_CGAL
138 }
139 // ----------------------------------------------------------------------------
140 // instanciation in library
141 // ----------------------------------------------------------------------------
142 #define _RHEOLEF_instanciation(T) \
143 template T orient2d ( \
144  const point_basic<T>&, \
145  const point_basic<T>&, \
146  const point_basic<T>&); \
147 template T orient3d ( \
148  const point_basic<T>&, \
149  const point_basic<T>&, \
150  const point_basic<T>&, \
151  const point_basic<T>&); \
152 template int sign_orient2d ( \
153  const point_basic<T>&, \
154  const point_basic<T>&, \
155  const point_basic<T>&); \
156 template int sign_orient3d ( \
157  const point_basic<T>&, \
158  const point_basic<T>&, \
159  const point_basic<T>&, \
160  const point_basic<T>&); \
161 
163 
164 } // namespace rheolef
see the Float page for the full documentation
point_basic< T >
Definition: piola_fem.h:135
rheolef::std value
Expr1::float_type T
Definition: field_expr.h:261
This file is part of Rheolef.
T orient2d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c)
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)
int sign_orient2d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c)
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
Float sgn(Float x)
Definition: sgn.icc:25
Float epsilon