|
Medusa
1.1
Coordinate Free Mehless Method implementation
|
|
Go to the documentation of this file. 1 #ifndef MEDUSA_BITS_DOMAINS_HALFLINKSREFINE_HPP_
2 #define MEDUSA_BITS_DOMAINS_HALFLINKSREFINE_HPP_
18 template <
typename vec_t>
24 template <
typename vec_t>
32 [&](
int i,
int j) { return domain.type(i) < domain.type(j); });
37 template <
class vec_t>
41 int region_size =
region.size();
42 int N = domain.
size();
45 assert_msg(region_size > 0,
"The region to refine is empty.");
48 int num_new_points = 0;
49 for (
int i = 0; i < region_size; i++) {
51 const vec_t pos = domain.
pos(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);
59 for (
int j = 1; j < n; ++j) {
61 vec_t candidate = 0.5 * (domain.
pos(c) + domain.
pos(s));
64 if (domain.
type(c) == 0 || domain.
type(s) == 0)
continue;
65 int candidate_type = std::max(domain.
type(c), domain.
type(s));
67 if (candidate_type < 0) {
69 if (normal_vector.squaredNorm() < 1e-15) {
73 normal_vector.normalize();
74 auto result = domain.
shape().projectPointToBoundary(candidate, normal_vector);
76 candidate = result.second;
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)
87 if (!domain.
shape().contains(candidate))
continue;
89 double dist2 = domain_tree.
query(candidate).second[0];
92 domain_tree.
insert(candidate);
93 if (candidate_type < 0) {
104 std::iota(new_ids.begin(), new_ids.end(), N);
110 #endif // MEDUSA_BITS_DOMAINS_HALFLINKSREFINE_HPP_
void insert(const vec_t &point)
Insert a point into the tree.
Root namespace for the whole library.
Range< int > region_
Node indexes around which to refine.
scalar_t dr(int i) const
Returns Euclidean distance to the second support node.
Scalar scalar_t
Type of the elements, alias of Scalar.
Class representing domain discretization along with an associated shape.
indexes_t all() const
Returns indexes of all nodes, i.e. {0, 1, ..., N-1}.
A k-d tree data structure that supports dynamic insertions and lazy-removal.
std::pair< mm::Range< int >, mm::Range< double > > query(const vec_t &point, int k=1)
Find k nearest neighbors to given point.
const Range< vec_t > & positions() const
Returns positions of all nodes.
#define assert_msg(cond,...)
Assert with better error reporting.
const vec_t & normal(int i) const
Returns outside unit normal of i-th node.
HalfLinksRefine & region(Range< int > region)
Set region to refine.
int size() const
Returns N, the number of nodes in this discretization.
static Range< int > refine_impl(DomainDiscretization< vec_t > &domain, const Range< int > ®ion, double fraction, KDTreeMutable< vec_t > &domain_tree)
Refine implementation with more control over fine grained details.
const DomainShape< vec_t > & shape() const
Returns geometric shape of the underlying domain.
int type(int i) const
Returns type of i-th node.
int addInternalNode(const vec_t &point, int type)
Adds a single interior node with specified type to this discretization.
const vec_t & pos(int i) const
Returns the position of i-th node.
Range< int > operator()(DomainDiscretization< vec_t > &domain) const
Refines given domain.
int size() const
Returns number of elements.
const Range< int > & support(int i) const
Returns support indices for i-th node.
container_t & sort(container_t &v)
Sorts a container inplace.
int addBoundaryNode(const vec_t &point, int type, const vec_t &normal)
Adds a boundary node with given type and normal to the domain.
HalfLinksRefine & min_dist(double fraction)
Minimal distance criterion around point p is that nodes generated at p must be at least fraction * cl...
double fraction_
Minimal distance fraction.