rheolef  6.6
geo_element_predicate.cc
Go to the documentation of this file.
1 #include "rheolef/geo_element.h"
2 
3 namespace rheolef {
4 
5 template <class T, class M>
6 bool
7 point_belongs_to_e (const array<point_basic<T>,M>& node, const geo_element& K, const point_basic<T>& x)
8 {
9  const point_basic<T>& p0 = node.dis_at(K[0]);
10  const point_basic<T>& p1 = node.dis_at(K[1]);
11  return x[0] >= p0[0] && x[0] <= p1[0];
12 }
13 template <class T, class M>
14 bool
15 point_belongs_to_t (const array<point_basic<T>,M>& node, const geo_element& K, const point_basic<T>& x)
16 {
17  static const T eps = 1e3*std::numeric_limits<T>::epsilon();
18  const point_basic<T>& p0 = node.dis_at(K[0]);
19  const point_basic<T>& p1 = node.dis_at(K[1]);
20  const point_basic<T>& p2 = node.dis_at(K[2]);
21  if (orient2d( x, p1, p2) < -eps) return false;
22  if (orient2d(p0, x, p2) < -eps) return false;
23  if (orient2d(p0, p1, x) < -eps) return false;
24  return true;
25 }
26 template <class T, class M>
27 bool
28 point_belongs_to_q (const array<point_basic<T>,M>& node, const geo_element& K, const point_basic<T>& x)
29 {
30  static const T eps = 1e3*std::numeric_limits<T>::epsilon();
31  const point_basic<T>& p0 = node.dis_at(K[0]);
32  const point_basic<T>& p1 = node.dis_at(K[1]);
33  const point_basic<T>& p2 = node.dis_at(K[2]);
34  const point_basic<T>& p3 = node.dis_at(K[3]);
35  if (orient2d(x, p0, p1) < -eps) return false;
36  if (orient2d(x, p1, p2) < -eps) return false;
37  if (orient2d(x, p2, p3) < -eps) return false;
38  if (orient2d(x, p3, p0) < -eps) return false;
39  return true;
40 }
41 template <class T, class M>
42 bool
43 point_belongs_to_T (const array<point_basic<T>,M>& node, const geo_element& K, const point_basic<T>& x)
44 {
45  static const T eps = 1e3*std::numeric_limits<T>::epsilon();
46  const point_basic<T>& p0 = node.dis_at(K[0]);
47  const point_basic<T>& p1 = node.dis_at(K[1]);
48  const point_basic<T>& p2 = node.dis_at(K[2]);
49  const point_basic<T>& p3 = node.dis_at(K[3]);
50  if (orient3d( x, p1, p2, p3) < -eps) return false;
51  if (orient3d(p0, x, p2, p3) < -eps) return false;
52  if (orient3d(p0, p1, x, p3) < -eps) return false;
53  if (orient3d(p0, p1, p2, x) < -eps) return false;
54  return true;
55 }
56 template <class T, class M>
57 bool
58 point_belongs_to_H (const array<point_basic<T>,M>& node, const geo_element& K, const point_basic<T>& x)
59 {
61  static const T eps = 1e3*std::numeric_limits<T>::epsilon();
62  const point_basic<T>& p0 = node.dis_at(K[0]);
63  const point_basic<T>& p1 = node.dis_at(K[1]);
64  const point_basic<T>& p2 = node.dis_at(K[2]);
65  const point_basic<T>& p3 = node.dis_at(K[3]);
66  const point_basic<T>& p4 = node.dis_at(K[4]);
67  const point_basic<T>& p5 = node.dis_at(K[5]);
68  const point_basic<T>& p6 = node.dis_at(K[6]);
69  const point_basic<T>& p7 = node.dis_at(K[7]);
70  const point_basic<T>* q[8] = {&p0, &p1, &p2, &p3, &p4, &p5, &p6, &p7};
71  for (size_type loc_isid = 0, loc_nsid = 6; loc_isid < loc_nsid; loc_isid++) {
72  size_type j0 = reference_element_H::subgeo_local_node (1, 2, loc_isid, 0);
73  size_type j1 = reference_element_H::subgeo_local_node (1, 2, loc_isid, 1);
74  size_type j2 = reference_element_H::subgeo_local_node (1, 2, loc_isid, 2);
75  if (orient3d(x, *(q[j0]), *(q[j1]), *(q[j2])) < -eps) return false;
76  }
77  return true;
78 }
79 template <class T, class M>
80 bool
81 point_belongs_to_P (const array<point_basic<T>,M>& node, const geo_element& K, const point_basic<T>& x)
82 {
84  static const T eps = 1e3*std::numeric_limits<T>::epsilon();
85  const point_basic<T>& p0 = node.dis_at(K[0]);
86  const point_basic<T>& p1 = node.dis_at(K[1]);
87  const point_basic<T>& p2 = node.dis_at(K[2]);
88  const point_basic<T>& p3 = node.dis_at(K[3]);
89  const point_basic<T>& p4 = node.dis_at(K[4]);
90  const point_basic<T>& p5 = node.dis_at(K[5]);
91  const point_basic<T>* q[6] = {&p0, &p1, &p2, &p3, &p4, &p5};
92  for (size_type loc_isid = 0, loc_nsid = 5; loc_isid < loc_nsid; loc_isid++) {
93  size_type j0 = reference_element_P::subgeo_local_node (1, 2, loc_isid, 0);
94  size_type j1 = reference_element_P::subgeo_local_node (1, 2, loc_isid, 1);
95  size_type j2 = reference_element_P::subgeo_local_node (1, 2, loc_isid, 2);
96  if (orient3d(x, *(q[j0]), *(q[j1]), *(q[j2])) < -eps) return false;
97  }
98  return true;
99 }
100 template <class T, class M>
101 bool
102 geo_element::contains (const array<point_basic<T>,M>& node, const point_basic<T>& x) const
103 {
104  switch (variant()) {
105  case reference_element::e: return point_belongs_to_e (node, *this, x);
106  case reference_element::t: return point_belongs_to_t (node, *this, x);
107  case reference_element::q: return point_belongs_to_q (node, *this, x);
108  case reference_element::T: return point_belongs_to_T (node, *this, x);
109  case reference_element::P: return point_belongs_to_P (node, *this, x);
110  case reference_element::H: return point_belongs_to_H (node, *this, x);
111  default: error_macro ("unsupported element type '" << name() << "'"); return false;
112  }
113 }
114 template bool geo_element::contains (const array<point_basic<Float>,sequential>&, const point_basic<Float>&) const;
115 #ifdef _RHEOLEF_HAVE_MPI
116 template bool geo_element::contains (const array<point_basic<Float>,distributed>&, const point_basic<Float>&) const;
117 #endif // _RHEOLEF_HAVE_MPI
118 
119 } // namespace rheolef
bool point_belongs_to_q(const array< point_basic< T >, M > &node, const geo_element &K, const point_basic< T > &x)
static const variant_type t
#define error_macro(message)
Definition: compiler.h:100
bool point_belongs_to_t(const array< point_basic< T >, M > &node, const geo_element &K, const point_basic< T > &x)
bool point_belongs_to_e(const array< point_basic< T >, M > &node, const geo_element &K, const point_basic< T > &x)
static const variant_type q
static const variant_type P
T orient2d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c)
point - vertex of a mesh
Definition: point.h:22
std::vector< int >::size_type size_type
irheostream, orheostream - large data streams
Definition: compiler.h:7
static const variant_type T
variant_type variant() const
static size_type subgeo_local_node(size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
static const variant_type e
bool contains(const array< point_basic< T >, M > &node, const point_basic< T > &x) const
static size_type subgeo_local_node(size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
static const variant_type H
bool point_belongs_to_P(const array< point_basic< T >, M > &node, const geo_element &K, const point_basic< T > &x)
bool point_belongs_to_T(const array< point_basic< T >, M > &node, const geo_element &K, const point_basic< T > &x)
bool point_belongs_to_H(const array< point_basic< T >, M > &node, const geo_element &K, const point_basic< T > &x)
T orient3d(const point_basic< T > &a, const point_basic< T > &b, const point_basic< T > &c, const point_basic< T > &d)
Float epsilon
geo_element - element of a mesh
size_type & node(size_type loc_inod)
reference_element_H::size_type size_type