1 #ifndef MEDUSA_BITS_DOMAINS_FINDBALANCEDSUPPORT_HPP_
2 #define MEDUSA_BITS_DOMAINS_FINDBALANCEDSUPPORT_HPP_
21 template <
class vec_t>
24 return normal.transpose().jacobiSvd(Eigen::ComputeFullV).matrixV().rightCols(
vec_t::dim-1);
28 template <
typename domain_t>
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()) {
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.",
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.",
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());
52 for (
int i : for_which) {
53 domain.support(i) = balancedSupport(
54 domain, tree, i, min_support_, max_support_, search_among, force_self_);
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) {
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) {
71 if (s[b] > 0) octant += (1 << b);
74 if (octant != -1 && !octants_covered[octant]) {
75 octants_covered[octant] =
true;
81 template <
typename domain_t,
typename kdtree_t>
83 int min_support_size,
int max_support_size,
84 const Range<int>& search_among,
bool force_self) {
86 typedef typename domain_t::vector_t vec_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));
92 int support_size = min_support_size;
93 const vec_t& center = domain.pos(i);
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);
101 while (pos_idx < max_support_size) {
102 indices = tree.query(center, support_size).first;
103 for (
int& j : indices) { j = search_among[j]; }
105 if (force_self && indices[0] != i) {
106 indices.insert(indices.begin(), i);
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--;
114 }
else if (domain.type(i) > 0) {
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);
123 if (support_size > max_support_size) { support_size = max_support_size; }
130 #endif // MEDUSA_BITS_DOMAINS_FINDBALANCEDSUPPORT_HPP_