You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 
 

125 lines
4.5 KiB

#ifndef MEDUSA_BITS_DOMAINS_GRAINDROPFILL_HPP_
#define MEDUSA_BITS_DOMAINS_GRAINDROPFILL_HPP_
/**
* @file
* Implementation of the GrainDropFill fill engine.
*/
#include "GrainDropFill_fwd.hpp"
#include <medusa/bits/utils/randutils.hpp>
#include <medusa/bits/spatial_search/KDTree.hpp>
#include <medusa/bits/spatial_search/Grid.hpp>
#include <medusa/bits/domains/discretization_helpers.hpp>
namespace mm {
template <typename vec_t>
GrainDropFill<vec_t>::GrainDropFill(const vec_t& bot, const vec_t& top) :
bot(bot), top(top), seed_(get_seed()) {
for (int i = 0; i < dim; ++i) {
assert_msg(top[i] > bot[i], "Bottom must be less than top, got bot = %s, top = %s.",
bot, top);
}
}
template <typename vec_t>
template <typename func_t>
void GrainDropFill<vec_t>::operator()(DomainDiscretization<vec_t>& domain, const func_t& h,
int type) const {
if (type == 0) type = 1;
std::vector<vec_t> nodes = fillBox(h);
KDTree<vec_t> boundary_tree(domain.positions());
for (const vec_t& p : nodes) {
if (!domain.contains(p, boundary_tree)) continue;
scalar_t r2 = boundary_tree.query(p).second[0];
scalar_t allowed_r = 0.75*h(p);
if (r2 >= allowed_r*allowed_r) {
domain.addInternalNode(p, type);
}
}
}
template <typename vec_t>
template <typename func_t>
std::vector<vec_t> GrainDropFill<vec_t>::fillBox(const func_t& h) const {
assert_msg(dx > 0, "Initial spacing must be positive, got %g. Did you set it correctly?", dx);
std::vector<vec_t> result;
// Initialization.
Vec<int, dim-1> grid_size;
for (int i = 0; i < dim-1; ++i) {
grid_size[i] = iceil((top[i] - bot[i]) / dx) + 1;
assert_msg(grid_size[i] >= 1, "Grid size in dimension %d is not positive.");
}
Grid<scalar_t, dim-1, int, Vec<int, dim-1>> z_coor(grid_size, bot[dim-1]);
Eigen::Map<Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>>
eig_z(z_coor.data().data(), z_coor.size());
eig_z += 1e-4*dx*Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>::Random(z_coor.size());
// Main iteration loop. Find locally minimal point and raise its neighbourhood up.
int min_lin_idx;
scalar_t min_z = eig_z.minCoeff(&min_lin_idx);
Vec<int, dim-1> min_idx = z_coor.multiIndex(min_lin_idx);
vec_t p;
int dot_nr = 0;
while (dot_nr < max_points && min_z <= top[dim-1] + excess_factor*dx) {
// Create new point.
for (int i = 0; i < dim-1; ++i) p[i] = bot[i] + min_idx[i] * dx;
p[dim-1] = min_z;
result.push_back(p);
++dot_nr;
// Update heights of neighbours.
scalar_t r = h(p);
int idx_rad = ifloor(r / dx);
Vec<int, dim-1> low = (min_idx - Vec<int, dim-1>::Constant(idx_rad)).cwiseMax(0);
Vec<int, dim-1> high = (min_idx + Vec<int, dim-1>::Constant(idx_rad+1)).cwiseMin(grid_size);
Vec<int, dim-1> cnt = low;
int cnt_lin;
do {
scalar_t height_update = r*r-((cnt-min_idx).template cast<scalar_t>()*dx).squaredNorm();
cnt_lin = z_coor.linearIndex(cnt);
if (height_update > 0) {
z_coor[cnt_lin] = std::max(z_coor[cnt_lin], min_z + std::sqrt(height_update));
}
} while (incrementCounter(cnt, low, high));
min_z = z_coor[min_lin_idx]; // Update minimum to be able to find something smaller later.
// Find next local minimum
int dist = 2*idx_rad;
while (dist > idx_rad) {
for (int i = 0; i < dim-1; ++i) {
low[i] = (((min_idx[i] - 2*idx_rad) % grid_size[i]) + grid_size[i]) % grid_size[i];
high[i] = (min_idx[i] + 2*idx_rad + 1) % grid_size[i];
}
auto old_min_idx = min_idx;
cnt = low;
do {
cnt_lin = z_coor.linearIndex(cnt);
if (z_coor[cnt_lin] < min_z) {
min_z = z_coor[cnt_lin];
min_lin_idx = cnt_lin;
min_idx = cnt;
}
} while (incrementCyclicCounter(cnt, low, high, grid_size));
dist = (min_idx - old_min_idx).cwiseAbs().maxCoeff();
}
}
return result;
}
/// Specialization for 1D
template <>
template <typename func_t>
std::vector<Vec1d> GrainDropFill<Vec1d>::fillBox(const func_t& h) const {
auto result = discretization_helpers::discretizeLineWithDensity(bot, top, h);
result.resize(std::min(result.size(), max_points));
return result;
}
} // namespace mm
#endif // MEDUSA_BITS_DOMAINS_GRAINDROPFILL_HPP_