#ifndef MEDUSA_BITS_DOMAINS_GENERALFILL_HPP_ #define MEDUSA_BITS_DOMAINS_GENERALFILL_HPP_ /** * @file * Implementation of the general node placing algorithm. */ #include "GeneralFill_fwd.hpp" #include "discretization_helpers.hpp" #include #include #include #include #include #include namespace mm { template GeneralFill::GeneralFill() : seed_(get_seed()) {} /// @cond template template void GeneralFill::operator()(DomainDiscretization &domain, const func_t &h, int type) const { KDTreeMutable tree; operator()(domain, h, tree, type); } template template void GeneralFill::operator()(DomainDiscretization &domain, const func_t &h, search_structure_t &search, int type) const { KDTree contains_search; domain.makeDiscreteContainsStructure(contains_search); auto contains_function = [&](const vec_t point) { return domain.contains(point, contains_search); }; operator()(domain, h, search, contains_function, type); } template template void GeneralFill::operator()(DomainDiscretization &domain, const func_t &h, search_structure_t &search, contains_function_t &contains_function, int type) const { if (type == 0) type = 1; std::mt19937 gen(seed_); std::cout << "Filling the domain interior...\n"; int cur_node = 0; int end_node = domain.size(); if (end_node == 0) { // If domain is empty, pick a random node inside it. vec_t lo_bound, hi_bound, random_node; std::tie(lo_bound, hi_bound) = domain.shape().bbox(); std::vector> distributions; for (int j = 0; j < dim; ++j) distributions.emplace_back(lo_bound[j], hi_bound[j]); int count = 0; do { for (int j = 0; j < dim; ++j) { random_node[j] = distributions[j](gen); } if (++count >= 10000) { std::string message = "No suitable node in domain could be found after 10000 tries." " This might happen if domain volume is very small compared " "to the volume of the bounding box. Try manually supplying " "the initial point."; throw std::runtime_error(message); } } while (!contains_function(random_node)); // std::cout << "新起始点:" << random_node[0] << ',' << random_node[1] << ',' << random_node[2] << '\n'; domain.addInternalNode(random_node, type); end_node = 1; } // Main algorithm loop. search.insert(domain.positions()); while ((cur_node < end_node) && (end_node < max_points)) { // std::cout << cur_node << " / " << end0_node << std::endl; vec_t p = domain.pos(cur_node); // std::cout << "新队头:" << p[0] << ',' << p[1] << ',' << p[2] << '\n'; scalar_t r = h(p); assert_msg(r > 0, "Nodal spacing radius should be > 0, got %g.", r); // std::cout << end_node << " : " << max_points << '\n'; auto candidates = discretization_helpers::SphereDiscretization::construct(r, n_samples, gen); // filter candidates regarding the domain and proximity criteria for (const auto &f : candidates) { vec_t node = p + f; // Shift center to point `p`. if (!contains_function(node)) continue; // If node is not in the domain. if (search.existsPointInSphere(node, r * zeta)) continue; // If node is not far enough // std::cout << "新候选点:" << node[0] << ',' << node[1] << ',' << node[2] << '\n'; domain.addInternalNode(node, type); // from other nodes. search.insert(node); end_node++; } cur_node++; } if (end_node >= max_points) { std::cerr << "Maximum number of points reached, fill may be incomplete." << std::endl; } std::cout << "Filling the domain interior done!\n"; } template GeneralFill &GeneralFill::proximityTolerance(scalar_t zeta) { assert_msg((0 < zeta && zeta < 1), "Zeta must be between 0 and 1, got %f.", zeta); this->zeta = zeta; return *this; } /// @endcond } // namespace mm #endif // MEDUSA_BITS_DOMAINS_GENERALFILL_HPP_