Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
point_predicate.cc
Go to the documentation of this file.
1
20// 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
37namespace rheolef {
38
39template<class T>
40static
41T
42inexact_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}
51template <class T>
52static
53T
54inexact_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}
70template <class T>
71int
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}
87template <class T>
88int
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}
105template<class T>
106T
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
119 return sgn*std::numeric_limits<T>::epsilon();
120#endif // _RHEOLEF_HAVE_CGAL
121}
122template <class T>
123T
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
136 return sgn*std::numeric_limits<T>::epsilon();
137#endif // _RHEOLEF_HAVE_CGAL
138}
139// ----------------------------------------------------------------------------
140// instanciation in library
141// ----------------------------------------------------------------------------
142#define _RHEOLEF_instanciation(T) \
143template T orient2d ( \
144 const point_basic<T>&, \
145 const point_basic<T>&, \
146 const point_basic<T>&); \
147template T orient3d ( \
148 const point_basic<T>&, \
149 const point_basic<T>&, \
150 const point_basic<T>&, \
151 const point_basic<T>&); \
152template int sign_orient2d ( \
153 const point_basic<T>&, \
154 const point_basic<T>&, \
155 const point_basic<T>&); \
156template 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
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
Expr1::float_type T
Definition field_expr.h:230
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_orient2d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c)
int sign_orient3d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c, const point_basic< T > &d)
Float sgn(Float x)
Definition sgn.icc:25