Medusa  1.1
Coordinate Free Mehless Method implementation
BasicRelax.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_BASICRELAX_HPP_
2 #define MEDUSA_BITS_DOMAINS_BASICRELAX_HPP_
3 
4 #include "BasicRelax_fwd.hpp"
5 #include <medusa/bits/utils/assert.hpp>
6 #include <medusa/bits/types/Range.hpp>
7 #include <medusa/bits/spatial_search/KDTree.hpp>
8 #include <medusa/bits/spatial_search/KDTreeMutable.hpp>
9 
15 namespace mm {
16 
17 
18 template <class domain_t, class radius_func_type>
19 void BasicRelax::operator()(domain_t& domain, const radius_func_type& r_func) const {
20  typedef typename domain_t::vector_t vec_t;
21  typedef typename domain_t::scalar_t scalar_t;
22  const int dim = domain_t::dim;
23 
24  // if no nodes are supplied, all interior nodes are processed
25  Range<int> nodes = (nodes_.empty()) ? static_cast<Range<int>>(domain.interior()) : nodes_;
26 
27  int N = domain.size();
28  assert_msg(N > num_neighbours, "Cannot relax a domain with less that num_neighbours nodes, "
29  "num_neighbours = %d, N = %d.", num_neighbours, N);
30  assert_msg(num_neighbours >= 1, "At least two neighbours are required, got %d.",
31  num_neighbours);
32  for (int i : nodes) {
33  assert_msg(0 <= i && i < N, "Node index %d out of range [0, %d) in relax.", i, N);
34  }
35 
36  Range<int> bnd = domain.boundary();
37  assert_msg(!bnd.empty(), "Relax requires boundary nodes, but this domain has none!");
38  KDTreeMutable<vec_t> boundary_tree(domain.positions()[bnd]);
39  KDTree<vec_t> tree;
40  Range<int> indices;
41 
42  Range<int> to_delete = {}; // nodes to be deleted from domain after relax
43  bool removed_nodes_warning = false;
44 
45  for (int c = 0; c < num_iterations; ++c) {
46  if (c == 0 || c % rebuild_tree_after == 0) {
47  tree.reset(domain.positions());
48  }
49  // nodes to be removed from relax procedure
50  Range<int> to_remove = {};
51  // Annealing factor -- depends only on iteration step
52  double a_f = ((initial_heat - final_heat) * (num_iterations - c) / num_iterations
53  + final_heat);
54  // Range of all displacements
55  Range<vec_t> dp(nodes.size(), vec_t(0.0));
56  int k;
57 // #if defined(_OPENMP)
58 // #pragma omp parallel for private(k) schedule(static)
59  // before enabling openmp vectors to_delete and to_remove
60  // have to be introduced for each thread and combined at the end of the loop
61  // otherwise expect se.g.\ fault
62 // #endif
63  for (k = 0; k < nodes.size(); ++k) {
64  int i = nodes[k];
65  assert_msg(domain.type(i) > 0, "Only interior nodes are to be relaxed");
66  vec_t pos = domain.pos(i);
67  indices = tree.query(pos, num_neighbours + 1).first;
68  double d0 = std::pow((domain.pos(indices[1]) - pos).norm(), 1.0 / dim);
69  // loop over support and compute potential contributions
70  scalar_t r_0 = 0;
71  for (int j = 1; j < indices.size(); ++j) {
72  // central node to support nodes vector
73  vec_t r = domain.pos(indices[j]) - pos;
74  assert_msg(r.norm() > 1e-15, "Nodes %d and %d have the same coordinates %s. "
75  "Try using lower initial heat.", i, indices[j], pos);
76  // potential gradient
77  // target density in terms of radius to the closest node
78  r_0 = r_func(domain.pos(indices[j]));
79  dp[k] += std::pow(std::abs(r_0) / r.norm(), potential_order) * r;
80  }
81  // no-distribution mode
82  if (r_0 < 0) dp[k] = 0.2 * dp[k].normalized() * a_f * d0;
83  dp[k] = -a_f * d0 * dp[k];
84  // Project escaped nodes on the boundary
85  if (!domain.shape().contains(pos + dp[k])) {
86  vec_t tmp, normal; // projection temp variables
87  bool project = true;
88  switch (projection_type) {
89  case DO_NOT_PROJECT: // mark escaped nodes for deletion
90  to_remove.push_back(k);
91  to_delete.push_back(nodes[k]);
92  project = false;
93  break;
94  case PROJECT_IN_DIRECTION: // project on boundary in direction of relax
95  tmp = 0.5 * (2 * pos + dp[k]);
96  normal = dp[k].normalized();
97  break;
98  case PROJECT_BETWEEN_CLOSEST: // project between closest nodes on boundary
99  Range<int> idx = boundary_tree.query(pos, 2).first;
100  int s = bnd[idx[0]], t = bnd[idx[1]];
101  tmp = 0.5 * (domain.pos(s) + domain.pos(t));
102  normal = (0.5 * (domain.normal(s) + domain.normal(t))).normalized();
103  break;
104  } // projection method selection end
105  if (project) {
106  // performs actual projection
107  bool success;
108  vec_t proj;
109  std::tie(success, proj) = domain.shape().projectPointToBoundary(tmp, normal);
110  if (success) {
111  // check if projected node is not too close to boundary node
112  Range<double> dist;
113  Range<int> ind;
114  std::tie(ind, dist) = boundary_tree.query(proj, 2);
115  if (std::sqrt(dist[0] / dist[1]) < boundary_projection_threshold) {
116  dp[k] = -dp[k];
117  continue;
118  }
119  double d_1 = std::sqrt(dist[0]);
120  double d_2 = std::sqrt(dist[1]);
121  double d_0 = std::sqrt(dist[0] + dist[1]);
122 
123  normal = ((1 - d_1 / d_0) * domain.normal(bnd[ind[0]])
124  + (1 - d_2 / d_0) * domain.normal(bnd[ind[1]])).normalized();
125  int boundary_type = domain.type(bnd[ind[0]]);
126  domain.changeToBoundary(i, proj, boundary_type, normal);
127  boundary_tree.insert(proj);
128  bnd.push_back(i);
129  // mark projected node to be removed from the list for relax
130  to_remove.push_back(k);
131  } else {
132  removed_nodes_warning = true;
133  // if projection failed, remove node
134  to_remove.push_back(k);
135  to_delete.push_back(nodes[k]);
136  }
137  } // project if end
138  } // escaped nodes if end
139  } // spatial loop end
140  nodes.remove(to_remove);
141  assert_msg(nodes.size() > num_neighbours + 1,
142  "No nodes in relax pool anymore, perhaps use lower heat");
143  // apply displacement
144  // #if defined(_OPENMP)
145  // #pragma omp parallel for private(k) schedule(static)
146  // #endif
147  for (k = 0; k < nodes.size(); ++k) {
148  domain.pos(nodes[k]) += dp[k];
149  }
150  // DEBUG OUT -- for plotting
151  // debug_out << "nodes(" << c + 1 << ",:,:)=" << domain.positions << ";\n";
152  }
153  // remove nodes at the end of relax to avoid mess with indices
154  domain.removeNodes(to_delete);
155  #if 0 // debug
156  debug_out << "nodes_final(:,:)=" << domain.positions << ";\n";
157  debug_out << "normal =" << domain.normals() << ";\n";
158  Range<int> boundary = domain.types < 0;
159  debug_out << "boundary =" << boundary << ";\n";
160  debug_out << "boundary_map =" << domain.boundary_map() << ";\n";
161  debug_out.close();
162  #endif
163  if (removed_nodes_warning) {
164  print_red("warning -- nodes removed during relax due to projection fail\n");
165  }
166 }
167 
168 } // namespace mm
169 
170 #endif // MEDUSA_BITS_DOMAINS_BASICRELAX_HPP_
mm::BasicRelax::nodes_
Range< int > nodes_
List of nodes to process.
Definition: BasicRelax_fwd.hpp:73
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
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::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::BasicRelax::projection_type
ProjectionType projection_type
On boundary projection method.
Definition: BasicRelax_fwd.hpp:74
dim
@ dim
Number of elements of this matrix.
Definition: MatrixBaseAddons.hpp:14
mm::BasicRelax::DO_NOT_PROJECT
@ DO_NOT_PROJECT
Escaped nodes are frozen outside the domain till the end of relax and then removed.
Definition: BasicRelax_fwd.hpp:51
mm::BasicRelax::PROJECT_IN_DIRECTION
@ PROJECT_IN_DIRECTION
Project on boundary in relax movement direction.
Definition: BasicRelax_fwd.hpp:56
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
mm::KDTree::query
std::pair< Range< int >, Range< scalar_t > > query(const vec_t &point, int k=1) const
Find k nearest neighbors to given point.
Definition: KDTree.hpp:16
mm::BasicRelax::operator()
void operator()(domain_t &domain) const
Runs the relax on the selected domain with constant distribution equals to domain characteristic dist...
Definition: BasicRelax_fwd.hpp:127
KDTreeMutable.hpp
mm::print_red
void print_red(const std::string &s)
Prints given text in bold red.
Definition: assert.hpp:85
mm::BasicRelax::boundary_projection_threshold
double boundary_projection_threshold
Threshold for projecting nodes on boundary.
Definition: BasicRelax_fwd.hpp:75
mm::BasicRelax::num_neighbours
int num_neighbours
Number of nodes to consider when calculating the potential.
Definition: BasicRelax_fwd.hpp:67
mm::BasicRelax::final_heat
double final_heat
Heat at the end of the relax, usually around 0.
Definition: BasicRelax_fwd.hpp:70
KDTree.hpp
mm::BasicRelax::initial_heat
double initial_heat
Initial heat, usually between 0 and 5.
Definition: BasicRelax_fwd.hpp:69
mm::BasicRelax::num_iterations
int num_iterations
Number of iterations performed.
Definition: BasicRelax_fwd.hpp:68
assert.hpp
mm::Range::size
int size() const
Returns number of elements.
Definition: Range_fwd.hpp:185
mm::BasicRelax::rebuild_tree_after
int rebuild_tree_after
How often engine rebuild search tree, 1 is perfect but slow.
Definition: BasicRelax_fwd.hpp:72
mm::BasicRelax::potential_order
int potential_order
Order of repulsing potential.
Definition: BasicRelax_fwd.hpp:71
BasicRelax_fwd.hpp
mm::BasicRelax::PROJECT_BETWEEN_CLOSEST
@ PROJECT_BETWEEN_CLOSEST
Project between two closest boundary nodes – this one might result in divergent behaviour,...
Definition: BasicRelax_fwd.hpp:63
mm::KDTree::reset
void reset(const Range< vec_t > &points)
Grows a new tree with new points.
Definition: KDTree_fwd.hpp:71
mm::Range< int >
mm::Range::remove
void remove(indexes_t indexes)
Remove elements wih given indexes.
Definition: Range.hpp:79
Range.hpp