Medusa  1.1
Coordinate Free Mehless Method implementation
GrainDropFill.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_GRAINDROPFILL_HPP_
2 #define MEDUSA_BITS_DOMAINS_GRAINDROPFILL_HPP_
3 
9 #include "GrainDropFill_fwd.hpp"
14 
15 namespace mm {
16 
17 template <typename vec_t>
18 GrainDropFill<vec_t>::GrainDropFill(const vec_t& bot, const vec_t& top) :
19  bot(bot), top(top), seed_(get_seed()) {
20  for (int i = 0; i < dim; ++i) {
21  assert_msg(top[i] > bot[i], "Bottom must be less than top, got bot = %s, top = %s.",
22  bot, top);
23  }
24 }
25 
26 template <typename vec_t>
27 template <typename func_t>
29  int type) const {
30  if (type == 0) type = 1;
31  std::vector<vec_t> nodes = fillBox(h);
32  KDTree<vec_t> boundary_tree(domain.positions());
33  for (const vec_t& p : nodes) {
34  if (!domain.contains(p, boundary_tree)) continue;
35  scalar_t r2 = boundary_tree.query(p).second[0];
36  scalar_t allowed_r = 0.75*h(p);
37  if (r2 >= allowed_r*allowed_r) {
38  domain.addInternalNode(p, type);
39  }
40  }
41 }
42 
43 template <typename vec_t>
44 template <typename func_t>
45 std::vector<vec_t> GrainDropFill<vec_t>::fillBox(const func_t& h) const {
46  assert_msg(dx > 0, "Initial spacing must be positive, got %g. Did you set it correctly?", dx);
47 
48  std::vector<vec_t> result;
49 
50  // Initialization.
51  Vec<int, dim-1> grid_size;
52  for (int i = 0; i < dim-1; ++i) {
53  grid_size[i] = iceil((top[i] - bot[i]) / dx) + 1;
54  assert_msg(grid_size[i] >= 1, "Grid size in dimension %d is not positive.");
55  }
56 
57  Grid<scalar_t, dim-1, int, Vec<int, dim-1>> z_coor(grid_size, bot[dim-1]);
58  Eigen::Map<Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>>
59  eig_z(z_coor.data().data(), z_coor.size());
60  eig_z += 1e-4*dx*Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>::Random(z_coor.size());
61 
62  // Main iteration loop. Find locally minimal point and raise its neighbourhood up.
63  int min_lin_idx;
64  scalar_t min_z = eig_z.minCoeff(&min_lin_idx);
65  Vec<int, dim-1> min_idx = z_coor.multiIndex(min_lin_idx);
66  vec_t p;
67  int dot_nr = 0;
68  while (dot_nr < max_points && min_z <= top[dim-1] + excess_factor*dx) {
69  // Create new point.
70  for (int i = 0; i < dim-1; ++i) p[i] = bot[i] + min_idx[i] * dx;
71  p[dim-1] = min_z;
72  result.push_back(p);
73  ++dot_nr;
74 
75  // Update heights of neighbours.
76  scalar_t r = h(p);
77  int idx_rad = ifloor(r / dx);
78  Vec<int, dim-1> low = (min_idx - Vec<int, dim-1>::Constant(idx_rad)).cwiseMax(0);
79  Vec<int, dim-1> high = (min_idx + Vec<int, dim-1>::Constant(idx_rad+1)).cwiseMin(grid_size);
80  Vec<int, dim-1> cnt = low;
81  int cnt_lin;
82  do {
83  scalar_t height_update = r*r-((cnt-min_idx).template cast<scalar_t>()*dx).squaredNorm();
84  cnt_lin = z_coor.linearIndex(cnt);
85  if (height_update > 0) {
86  z_coor[cnt_lin] = std::max(z_coor[cnt_lin], min_z + std::sqrt(height_update));
87  }
88  } while (incrementCounter(cnt, low, high));
89  min_z = z_coor[min_lin_idx]; // Update minimum to be able to find something smaller later.
90 
91  // Find next local minimum
92  int dist = 2*idx_rad;
93  while (dist > idx_rad) {
94  for (int i = 0; i < dim-1; ++i) {
95  low[i] = (((min_idx[i] - 2*idx_rad) % grid_size[i]) + grid_size[i]) % grid_size[i];
96  high[i] = (min_idx[i] + 2*idx_rad + 1) % grid_size[i];
97  }
98  auto old_min_idx = min_idx;
99  cnt = low;
100  do {
101  cnt_lin = z_coor.linearIndex(cnt);
102  if (z_coor[cnt_lin] < min_z) {
103  min_z = z_coor[cnt_lin];
104  min_lin_idx = cnt_lin;
105  min_idx = cnt;
106  }
107  } while (incrementCyclicCounter(cnt, low, high, grid_size));
108  dist = (min_idx - old_min_idx).cwiseAbs().maxCoeff();
109  }
110  }
111  return result;
112 }
113 
115 template <>
116 template <typename func_t>
117 std::vector<Vec1d> GrainDropFill<Vec1d>::fillBox(const func_t& h) const {
118  auto result = discretization_helpers::discretizeLineWithDensity(bot, top, h);
119  result.resize(std::min(result.size(), max_points));
120  return result;
121 }
122 
123 } // namespace mm
124 
125 #endif // MEDUSA_BITS_DOMAINS_GRAINDROPFILL_HPP_
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::incrementCyclicCounter
bool incrementCyclicCounter(Vec< int, dim > &counter, const Vec< int, dim > &low, const Vec< int, dim > &high, const Vec< int, dim > &size)
Increments a multi-dimensional counter with given upper and lower limits and global upper size,...
Definition: numutils.hpp:182
mm::DomainDiscretization
Class representing domain discretization along with an associated shape.
Definition: DomainDiscretization_fwd.hpp:46
mm::Vec
Eigen::Matrix< scalar_t, dim, 1, Eigen::ColMajor|Eigen::AutoAlign, dim, 1 > Vec
Fixed size vector type, representing a mathematical 1d/2d/3d vector.
Definition: Vec_fwd.hpp:31
mm::KDTree
Class representing a static k-d tree data structure.
Definition: KDTree_fwd.hpp:36
mm::GrainDropFill::bot
vec_t bot
Bottom left corner of the box to fill.
Definition: GrainDropFill_fwd.hpp:46
dim
@ dim
Number of elements of this matrix.
Definition: MatrixBaseAddons.hpp:14
mm::GrainDropFill::fillBox
std::vector< vec_t > fillBox(const func_t &h) const
Actual fill algorithm which fills [bot, top] box with spacing h.
Definition: GrainDropFill.hpp:45
randutils.hpp
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::GrainDropFill::dim
@ dim
Dimensionality of the domain.
Definition: GrainDropFill_fwd.hpp:42
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::ifloor
int ifloor(T x)
Floors a floating point to an integer.
Definition: numutils.hpp:36
mm::GrainDropFill::scalar_t
vec_t::scalar_t scalar_t
Scalar type;.
Definition: GrainDropFill_fwd.hpp:39
mm::iceil
int iceil(T x)
Ceils a floating point to an integer.
Definition: numutils.hpp:27
mm::Grid
Class representing a simple n-dimensional grid structure, which supports indexing and storing values.
Definition: Grid_fwd.hpp:36
mm::GrainDropFill::operator()
void operator()(DomainDiscretization< vec_t > &domain, const func_t &h, int type=0) const
Fills domain with a quality node distribution.
Definition: GrainDropFill.hpp:28
mm::DomainDiscretization::contains
bool contains(const vec_t &point) const
Returns true if point is inside the domain.
Definition: DomainDiscretization.hpp:315
mm::get_seed
unsigned int get_seed()
Return a random seed.
Definition: randutils.hpp:25
KDTree.hpp
mm::discretization_helpers::discretizeLineWithDensity
Range< vec_t > discretizeLineWithDensity(const vec_t &p, const vec_t &q, const func_t &delta_r)
Returns nodes lying on the line segment pq with approximate distances delta_r.
Definition: discretization_helpers.hpp:35
Grid.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::GrainDropFill::top
vec_t top
Top right corner of the box to fill.
Definition: GrainDropFill_fwd.hpp:47
discretization_helpers.hpp
GrainDropFill_fwd.hpp
mm::incrementCounter
bool incrementCounter(Vec< int, dim > &counter, const Vec< int, dim > &limit)
Increments a multi-dimensional counter with given limits.
Definition: numutils.hpp:120
mm::GrainDropFill::GrainDropFill
GrainDropFill(const vec_t &bot, const vec_t &top)
Prepare to fill domain within [bot, top] bounding box.
Definition: GrainDropFill.hpp:18