Medusa  1.1
Coordinate Free Mehless Method implementation
HalfLinksRefine.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_HALFLINKSREFINE_HPP_
2 #define MEDUSA_BITS_DOMAINS_HALFLINKSREFINE_HPP_
3 
10 #include "HalfLinksRefine_fwd.hpp"
14 #include <numeric>
15 
16 namespace mm {
17 
18 template <typename vec_t>
20  KDTreeMutable<vec_t> domain_tree(domain.positions());
21  return operator()(domain, domain_tree);
22 }
23 
24 template <typename vec_t>
28  if (region.empty()) { region = domain.all(); }
29 
30  // sort: boundary nodes first
31  std::sort(region.begin(), region.end(),
32  [&](int i, int j) { return domain.type(i) < domain.type(j); });
33 
34  return refine_impl(domain, region, fraction_, tree);
35 }
36 
37 template <class vec_t>
39  DomainDiscretization<vec_t>& domain, const Range<int>& region,
40  double fraction, KDTreeMutable<vec_t>& domain_tree) {
41  int region_size = region.size();
42  int N = domain.size();
43  using scalar_t = typename vec_t::scalar_t;
44 
45  assert_msg(region_size > 0, "The region to refine is empty.");
46 
47  // Iterate through points in region and generate new points
48  int num_new_points = 0;
49  for (int i = 0; i < region_size; i++) {
50  int c = region[i]; // the global domain index
51  const vec_t pos = domain.pos(c);
52  const Range<int> supp = domain.support(c);
53  assert_msg(supp.size() >= 2, "At least 2 nodes must be in support of every node, %d "
54  "found in support of node %d.", supp.size(), c);
55  scalar_t min_dist = fraction * domain.dr(c);
56 
57  int n = supp.size();
58  // Half links to my support.
59  for (int j = 1; j < n; ++j) {
60  int s = supp[j];
61  vec_t candidate = 0.5 * (domain.pos(c) + domain.pos(s));
62 
63  // Decide the type of new node and project to boundary if necessary.
64  if (domain.type(c) == 0 || domain.type(s) == 0) continue;
65  int candidate_type = std::max(domain.type(c), domain.type(s));
66  vec_t normal_vector;
67  if (candidate_type < 0) {
68  normal_vector = domain.normal(c) + domain.normal(s);
69  if (normal_vector.squaredNorm() < 1e-15) {
70  // normal_vector of given points point in opposite directions.
71  continue;
72  }
73  normal_vector.normalize();
74  auto result = domain.shape().projectPointToBoundary(candidate, normal_vector);
75  if (result.first) {
76  candidate = result.second;
77  } else {
78  std::cerr << format("Adding point %s with type %d to the boundary along "
79  "normal %s was not successful.",
80  candidate, candidate_type, normal_vector)
81  << std::endl;
82  continue;
83  }
84  }
85 
86  // If nodes are out of the domain, they should be thrown away.
87  if (!domain.shape().contains(candidate)) continue;
88  // Distance to closest node.
89  double dist2 = domain_tree.query(candidate).second[0];
90  if (dist2 >= min_dist * min_dist) { // add new point
91  ++num_new_points;
92  domain_tree.insert(candidate);
93  if (candidate_type < 0) {
94  domain.addBoundaryNode(candidate, candidate_type, normal_vector);
95  } else {
96  domain.addInternalNode(candidate, candidate_type);
97  }
98  }
99  }
100  }
101 
102  // return back indices of new nodes
103  Range<int> new_ids(num_new_points);
104  std::iota(new_ids.begin(), new_ids.end(), N);
105  return new_ids;
106 }
107 
108 } // namespace mm
109 
110 #endif // MEDUSA_BITS_DOMAINS_HALFLINKSREFINE_HPP_
mm::KDTreeMutable::insert
void insert(const vec_t &point)
Insert a point into the tree.
Definition: KDTreeMutable.hpp:24
mm
Root namespace for the whole library.
Definition: Gaussian.hpp:14
mm::HalfLinksRefine::region_
Range< int > region_
Node indexes around which to refine.
Definition: HalfLinksRefine_fwd.hpp:30
mm::DomainDiscretization::dr
scalar_t dr(int i) const
Returns Euclidean distance to the second support node.
Definition: DomainDiscretization_fwd.hpp:148
scalar_t
Scalar scalar_t
Type of the elements, alias of Scalar.
Definition: MatrixBaseAddons.hpp:16
mm::DomainDiscretization
Class representing domain discretization along with an associated shape.
Definition: DomainDiscretization_fwd.hpp:46
mm::DomainDiscretization::all
indexes_t all() const
Returns indexes of all nodes, i.e. {0, 1, ..., N-1}.
Definition: DomainDiscretization_fwd.hpp:186
mm::KDTreeMutable
A k-d tree data structure that supports dynamic insertions and lazy-removal.
Definition: HalfLinksRefine_fwd.hpp:18
mm::KDTreeMutable::query
std::pair< mm::Range< int >, mm::Range< double > > query(const vec_t &point, int k=1)
Find k nearest neighbors to given point.
Definition: KDTreeMutable.hpp:45
mm::DomainDiscretization::positions
const Range< vec_t > & positions() const
Returns positions of all nodes.
Definition: DomainDiscretization_fwd.hpp:127
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
mm::DomainDiscretization::normal
const vec_t & normal(int i) const
Returns outside unit normal of i-th node.
Definition: DomainDiscretization.hpp:40
mm::HalfLinksRefine::region
HalfLinksRefine & region(Range< int > region)
Set region to refine.
Definition: HalfLinksRefine_fwd.hpp:38
mm::DomainDiscretization::size
int size() const
Returns N, the number of nodes in this discretization.
Definition: DomainDiscretization_fwd.hpp:189
KDTreeMutable.hpp
mm::HalfLinksRefine::refine_impl
static Range< int > refine_impl(DomainDiscretization< vec_t > &domain, const Range< int > &region, double fraction, KDTreeMutable< vec_t > &domain_tree)
Refine implementation with more control over fine grained details.
Definition: HalfLinksRefine.hpp:38
mm::DomainDiscretization::shape
const DomainShape< vec_t > & shape() const
Returns geometric shape of the underlying domain.
Definition: DomainDiscretization_fwd.hpp:162
mm::DomainDiscretization::type
int type(int i) const
Returns type of i-th node.
Definition: DomainDiscretization_fwd.hpp:158
HalfLinksRefine_fwd.hpp
assert.hpp
mm::DomainDiscretization::addInternalNode
int addInternalNode(const vec_t &point, int type)
Adds a single interior node with specified type to this discretization.
Definition: DomainDiscretization.hpp:49
mm::DomainDiscretization::pos
const vec_t & pos(int i) const
Returns the position of i-th node.
Definition: DomainDiscretization_fwd.hpp:129
mm::HalfLinksRefine::operator()
Range< int > operator()(DomainDiscretization< vec_t > &domain) const
Refines given domain.
Definition: HalfLinksRefine.hpp:19
mm::Range::size
int size() const
Returns number of elements.
Definition: Range_fwd.hpp:185
mm::DomainDiscretization::support
const Range< int > & support(int i) const
Returns support indices for i-th node.
Definition: DomainDiscretization_fwd.hpp:138
mm::sort
container_t & sort(container_t &v)
Sorts a container inplace.
Definition: stdtypesutils.hpp:43
mm::Range< int >
Range.hpp
mm::DomainDiscretization::addBoundaryNode
int addBoundaryNode(const vec_t &point, int type, const vec_t &normal)
Adds a boundary node with given type and normal to the domain.
Definition: DomainDiscretization.hpp:56
mm::HalfLinksRefine::min_dist
HalfLinksRefine & min_dist(double fraction)
Minimal distance criterion around point p is that nodes generated at p must be at least fraction * cl...
Definition: HalfLinksRefine_fwd.hpp:44
mm::HalfLinksRefine::fraction_
double fraction_
Minimal distance fraction.
Definition: HalfLinksRefine_fwd.hpp:31