Medusa  1.1
Coordinate Free Mehless Method implementation
FindBalancedSupport.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_FINDBALANCEDSUPPORT_HPP_
2 #define MEDUSA_BITS_DOMAINS_FINDBALANCEDSUPPORT_HPP_
3 
11 #include "Eigen/SVD"
12 
18 namespace mm {
19 
21 template <class vec_t>
23 FindBalancedSupport::getFrame(const vec_t& normal) {
24  return normal.transpose().jacobiSvd(Eigen::ComputeFullV).matrixV().rightCols(vec_t::dim-1);
25 }
27 
28 template <typename domain_t>
29 void mm::FindBalancedSupport::operator()(domain_t& domain) const {
30  auto for_which = for_which_;
31  if (for_which.empty()) for_which = domain.all();
32  auto search_among = search_among_;
33  if (search_among.empty()) {
34  search_among = Range<int>::seq(domain.size()); // include potential ghost nodes
35  }
36  assert_msg(!domain.positions().empty(), "Cannot find support in an empty domain.");
37  assert_msg(min_support_ > 0, "Support size must be greater than 0, got %d.",
38  min_support_);
39  assert_msg(min_support_ <= search_among.size(), "Support size (%d) cannot exceed number of "
40  "points that we are searching among (%d).", min_support_, search_among.size());
41  assert_msg(!for_which.empty(), "Set of nodes for which to find the support is empty.");
42  assert_msg(!search_among.empty(), "Set of nodes to search among is empty.");
43  for (int x : for_which) {
44  assert_msg(0 <= x && x < domain.size(), "Index %d out of range [%d, %d) in for_which.",
45  x, 0, domain.size());
46  }
47  for (int x : search_among) {
48  assert_msg(0 <= x && x < domain.size(), "Index %d out of range [%d, %d) in "
49  "search_among.", x, 0, domain.size());
50  }
51  KDTree<typename domain_t::vector_t> tree(domain.positions()[search_among]);
52  for (int i : for_which) {
53  domain.support(i) = balancedSupport(
54  domain, tree, i, min_support_, max_support_, search_among, force_self_);
55  }
56 }
57 
58 template <class vec_t>
60  const vec_t& dx, typename vec_t::scalar_t tol, std::vector<bool>& octants_covered,
61  const Eigen::Matrix<typename vec_t::scalar_t, vec_t::dim, -1>& basis) {
62  // for each basis vector, determine if it is in the same or opposite direction
63  // record that value in `octants`
64  int octant = 0;
65  Eigen::Matrix<typename vec_t::scalar_t, Eigen::Dynamic, 1> s = basis.transpose()*dx;
66  for (int b = 0; b < basis.cols(); ++b) {
67  if (std::abs(s[b]) < tol) {
68  octant = -1;
69  break;
70  }
71  if (s[b] > 0) octant += (1 << b);
72  }
73  // mark computed octant as covered
74  if (octant != -1 && !octants_covered[octant]) {
75  octants_covered[octant] = true;
76  return true;
77  }
78  return false;
79 }
80 
81 template <typename domain_t, typename kdtree_t>
82 Range<int> FindBalancedSupport::balancedSupport(domain_t& domain, const kdtree_t& tree, int i,
83  int min_support_size, int max_support_size,
84  const Range<int>& search_among, bool force_self) {
85  const int dim = domain_t::dim;
86  typedef typename domain_t::vector_t vec_t;
87  typedef typename domain_t::scalar_t scalar_t;
88  Eigen::Matrix<scalar_t, dim, Eigen::Dynamic> basis =
89  Eigen::Matrix<scalar_t, dim, dim>::Identity();
90  if (domain.type(i) < 0) basis = getFrame(domain.normal(i));
91 
92  int support_size = min_support_size;
93  const vec_t& center = domain.pos(i);
94  double tol = 1e-6;
95  int target = 1 << basis.cols();
96  int boundary_neighbours_required = (domain.type(i) < 0) ? 2*(dim-1)+1 : 0;
97  std::vector<bool> octants_covered(target, false);
98  int pos_idx = 0;
99  Range<int> indices;
100 
101  while (pos_idx < max_support_size) {
102  indices = tree.query(center, support_size).first;
103  for (int& j : indices) { j = search_among[j]; }
104 
105  if (force_self && indices[0] != i) {
106  indices.insert(indices.begin(), i);
107  }
108 
109  for (; pos_idx < indices.size(); ++pos_idx) {
110  vec_t dx = domain.pos(indices[pos_idx]) - center;
111  if (domain.type(indices[pos_idx]) < 0 && domain.type(i) < 0) {
112  boundary_neighbours_required--;
113  target -= mark_quadrant(dx, tol, octants_covered, basis);
114  } else if (domain.type(i) > 0) {
115  target -= mark_quadrant(dx, tol, octants_covered, basis);
116  }
117  if (target <= 0 && boundary_neighbours_required <= 0 &&
118  pos_idx >= min_support_size-1) {
119  return Range<int>(indices.begin(), indices.begin()+pos_idx+1);
120  }
121  }
122  support_size *= 2;
123  if (support_size > max_support_size) { support_size = max_support_size; }
124  }
125  return indices;
126 }
127 
128 } // namespace mm
129 
130 #endif // MEDUSA_BITS_DOMAINS_FINDBALANCEDSUPPORT_HPP_
mm::FindBalancedSupport::balancedSupport
static Range< int > balancedSupport(domain_t &domain, const kdtree_t &tree, int i, int min_support_size, int max_support_size, const Range< int > &search_among, bool force_self)
Finds support for a given node i.
Definition: FindBalancedSupport.hpp:82
mm
Root namespace for the whole library.
Definition: Gaussian.hpp:14
scalar_t
Scalar scalar_t
Type of the elements, alias of Scalar.
Definition: MatrixBaseAddons.hpp:16
mm::KDTree
Class representing a static k-d tree data structure.
Definition: KDTree_fwd.hpp:36
mm::FindBalancedSupport::getFrame
static Eigen::Matrix< typename vec_t::scalar_t, vec_t::dim, vec_t::dim-1 > getFrame(const vec_t &normal)
Gets orthonormal basis of vectors perpendicular to normal.
dim
@ dim
Number of elements of this matrix.
Definition: MatrixBaseAddons.hpp:14
mm::FindBalancedSupport::operator()
void operator()(domain_t &domain) const
Find support for nodes in domain.
Definition: FindBalancedSupport.hpp:29
mm::FindBalancedSupport::mark_quadrant
static bool mark_quadrant(const vec_t &dx, typename vec_t::scalar_t tol, std::vector< bool > &octants_covered, const Eigen::Matrix< typename vec_t::scalar_t, vec_t::dim, -1 > &basis)
Marks a quadrant in which dx resides as occupied.
Definition: FindBalancedSupport.hpp:59
mm::Range::seq
static Range seq(V n)
Returns range {0, ..., n-1}.
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
FindBalancedSupport_fwd.hpp
KDTree.hpp
Matrix
Matrix(const Scalar &s)
Construct matrix from scalar. Enabled only for fixed size matrices.
Definition: MatrixAddons.hpp:21
mm::Range< int >