Homec4science
Diffusion CADD_mesher (dev/incorporate-qhull) Browse

rCADDMESH/src/common/point_container.cc
58a104d7bd95dev/incorporate-qhull

Last Change
Raw File

point_container.cc
View Options

#include "point_container.hh"
#include <cmath>
#include <vector>
#include <algorithm>
#include "geometry.hh"
#include <stdexcept>
#include <set>
//template <Uint DIM>
//using IntCoord = std::array<int, DIM>;
//template <Uint DIM>
//using GridType = std::map<std::array<int, DIM>, std::vector<Uint>>;
//template <Uint DIM>
//using RevGridType = std::map<Uint, std::array<int, DIM>>;
template<Uint DIM>
void PointContainer<DIM>::filter(const Geometry<DIM> & geom,
Real tol) {
PointContainer temp;
for (Uint point = 0 ; point < this->getSize() ; ++point) {
PointRef<DIM> current_point = this->getPoint(point);
if (geom.is_inside(current_point, tol)) {
temp.addPoint(current_point);
}
}
this->values = std::move(temp.values);
this->modified = true;
}
// accumulates all the points of my box and all neighbouring boxes
template<Uint DIM>
inline void accumulate_test_points(const Uint & direction,
const std::array<int, DIM> & key,
std::map<std::array<int, DIM>,
std::vector<Uint>> & grid,
std::array<int, DIM> & off_set,
std::vector<Uint> & test_points) {
for (int i = -1 ; i < 2 ; ++i) {
off_set[direction] = key[direction] + i;
if (direction < DIM -1) {
accumulate_test_points<DIM>(direction + 1, key, grid, off_set, test_points);
} else {
for (Uint point_id: grid[off_set]) {
test_points.push_back(point_id);
}
}
}
}
template<Uint DIM>
void clean_grid(const Uint & index,
std::map<std::array<int, DIM>, std::vector<Uint>> & grid,
std::map<Uint, std::array<int, DIM>> & reverse_grid) {
std::array<int, DIM> & box = reverse_grid[index];
auto pos = std::find(grid.at(box).begin(), grid.at(box).end(), index);
grid.at(box).erase(pos);
reverse_grid.erase(index);
}
template<Uint DIM>
void PointContainer<DIM>::clear() {
this->threshold = -1;
this->values.clear();
this->modified = true;
}
template<Uint DIM>
Uint PointContainer<DIM>::reportDuplicates(Real threshold) {
this->setThreshold(threshold);
return this->sniffDuplicates<false>();
}
template<Uint DIM>
Uint PointContainer<DIM>::cleanDuplicates(Real threshold) {
this->setThreshold(threshold);
return this->sniffDuplicates<true>();
}
template<Uint DIM>
void PointContainer<DIM>::createSearchGrid() {
this->grid.clear();
this->reverse_grid.clear();
this->boxes.clear();
// sort all the points in a grid of size threshold
std::array<int, DIM> tmp_point;
for (Uint point_id = 0 ; point_id < this->getSize() ; ++point_id) {
for (Uint dim = 0 ; dim < DIM ; ++dim) {
tmp_point[dim] = std::floor(this->getConstPoint(point_id)[dim]/this->threshold);
}
this->grid[tmp_point].push_back(point_id);
this->reverse_grid[point_id] = tmp_point;
}
for (auto item: grid) {
this->boxes.emplace_back(item.first);
}
}
template<Uint DIM>
bool PointContainer<DIM>::containsExacly(const PointRef<DIM> & testpoint){
if (this->threshold <= 0) {
this->setThreshold(1);
}
std::array<int, DIM> tmp_point;
for (Uint dim = 0 ; dim < DIM ; ++dim) {
tmp_point[dim] = std::floor(testpoint[dim]/this->threshold);
}
for (auto point_id: this->grid[tmp_point]) {
if (this->getConstPoint(point_id) == testpoint) {
return true;
}
}
return false;
}
template<Uint DIM>
void PointContainer<DIM>::setThreshold(Real thresh)
{
if (thresh <= 0.) {
throw std::runtime_error("The threshold needs to be a positive finite real scalar");
}
if ((thresh != this->threshold) || this->modified) {
this->threshold = thresh;
this->createSearchGrid();
}
}
template<Uint DIM>
template<bool do_delete>
Uint PointContainer<DIM>::sniffDuplicates(){
Real tol2 = this->threshold*this->threshold;
std::array<int, DIM> tmp_point;
// get potential neighbours for each point in the container
std::set<Uint> nodes_to_delete;
Uint counter = 0;
for (auto key: this->boxes) {
std::vector<Uint> test_points;
accumulate_test_points<DIM>(0, key, this->grid, tmp_point, test_points);
std::vector<Uint> & local_ids = this->grid[key];
for (Uint me : local_ids) {
if (nodes_to_delete.find(me) == nodes_to_delete.end()) {
for (Uint other:test_points) {
if ((me != other) &&
(nodes_to_delete.find(other) == nodes_to_delete.end())) {
Real r_square = distance2(this->getConstPoint(me),
this->getConstPoint(other));
if (r_square < tol2) {
//There's a violation! kill it
counter++;
clean_grid<DIM>(other, this->grid, this->reverse_grid);
if (do_delete) {
std::cout << "Deleting point " << other << " at position" << this->getConstPoint(other) << std::endl;
nodes_to_delete.insert(other);
} else {
std::cout << "Duplicates: " << other << " at ";
this->getConstPoint(other).print_coords(std::cout);
std::cout << std::endl;
std::cout << " and: " << me << " at ";
this->getConstPoint(me).print_coords(std::cout);
std::cout << std::endl;
std::cout << " distance: " << std::sqrt(r_square) << std::endl;
}
}
}
}
}
}
}
if (do_delete) {
if (counter > 0 ) {
this->modified = true;
for (auto it = nodes_to_delete.rbegin(); it != nodes_to_delete.rend();
++it) {
this->deletePoint(*it);
}
}
}
return counter;
}
template class PointContainer<twoD>;
template class PointContainer<threeD>;
c4science · Help