← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to nfem/geo_element/point_predicate.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/geo_element/point_predicate.cc
 
1
/// This file is part of Rheolef.
 
2
///
 
3
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
4
///
 
5
/// Rheolef is free software; you can redistribute it and/or modify
 
6
/// it under the terms of the GNU General Public License as published by
 
7
/// the Free Software Foundation; either version 2 of the License, or
 
8
/// (at your option) any later version.
 
9
///
 
10
/// Rheolef is distributed in the hope that it will be useful,
 
11
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
/// GNU General Public License for more details.
 
14
///
 
15
/// You should have received a copy of the GNU General Public License
 
16
/// along with Rheolef; if not, write to the Free Software
 
17
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
18
/// 
 
19
/// =========================================================================
 
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
 
 
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
 
72
sign_orient2d (
 
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
 
89
sign_orient3d (
 
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
static
 
107
T
 
108
orient2d (const point_basic<T>& a, const point_basic<T>& b,
 
109
          const point_basic<T>& c)
 
110
{
 
111
  T value = inexact_orient2d(a, b, c);
 
112
#ifndef _RHEOLEF_HAVE_CGAL
 
113
  return value;
 
114
#else // _RHEOLEF_HAVE_CGAL
 
115
  int sgn =    sign_orient2d(a, b, c);
 
116
  value = fabs(value);
 
117
  if (sgn == 0)   return 0;
 
118
  if (value != 0) return sgn*value;
 
119
  // sgn != 0 but value == 0
 
120
  return sgn*std::numeric_limits<T>::epsilon();
 
121
#endif // _RHEOLEF_HAVE_CGAL
 
122
}
 
123
template <class T>
 
124
static
 
125
T
 
126
orient3d (const point_basic<T>& a, const point_basic<T>& b,
 
127
          const point_basic<T>& c, const point_basic<T>& d)
 
128
{
 
129
  T value = inexact_orient3d(a, b, c, d);
 
130
#ifndef _RHEOLEF_HAVE_CGAL
 
131
  return value;
 
132
#else // _RHEOLEF_HAVE_CGAL
 
133
  int sgn =    sign_orient3d(a, b, c, d);
 
134
  value = fabs(value);
 
135
  if (sgn == 0)   return 0;
 
136
  if (value != 0) return sgn*value;
 
137
  // sgn != 0 but value == 0
 
138
  return sgn*std::numeric_limits<T>::epsilon();
 
139
#endif // _RHEOLEF_HAVE_CGAL
 
140
}
 
141
// ----------------------------------------------------------------------------
 
142
// instanciation in library
 
143
// ----------------------------------------------------------------------------
 
144
#define _RHEOLEF_instanciation(T)                                       \
 
145
template T orient2d (                                                   \
 
146
  const point_basic<T>&,                                                \
 
147
  const point_basic<T>&,                                                \
 
148
  const point_basic<T>&);                                               \
 
149
template T orient3d (                                                   \
 
150
  const point_basic<T>&,                                                \
 
151
  const point_basic<T>&,                                                \
 
152
  const point_basic<T>&,                                                \
 
153
  const point_basic<T>&);                                               \
 
154
template int sign_orient2d (                                            \
 
155
  const point_basic<T>&,                                                \
 
156
  const point_basic<T>&,                                                \
 
157
  const point_basic<T>&);                                               \
 
158
template int sign_orient3d (                                            \
 
159
  const point_basic<T>&,                                                \
 
160
  const point_basic<T>&,                                                \
 
161
  const point_basic<T>&,                                                \
 
162
  const point_basic<T>&);                                               \
 
163
 
 
164
_RHEOLEF_instanciation(Float)
 
165
 
 
166
} // namespace rheolef