rheolef  6.5
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
120