Homec4science
Diffusion CADD_mesher (dev/incorporate-qhull) Browse

rCADDMESH/src/domain/tree.cc
cfac30bd70a4dev/incorporate-qhull

Last Change
Raw File

tree.cc
View Options

#include "tree.hh"
#include "common.hh"
#include <limits>
template <Uint DIM>
Tree<DIM>::Tree(const Real * centre_, Uint size_, Lattice<DIM> * lattice_, Domain<DIM> * domain_)
:size(size_), domain(domain_), lattice(lattice_), neighbours(0), level(0) {
this->nb_boxes ++;
if (centre_ == NULL) {
for (Uint i = 0; i < DIM ; ++i) {
this->south_west[i] = 0;
}
} else {
vector_set<DIM>(this->south_west, centre_);
}
for (Uint latt_dir = 0; latt_dir < DIM ; ++latt_dir) {
vector_incr<DIM>(this->south_west, this->lattice->getLatticeVector(latt_dir),
-0.5*this->size * this->lattice->getConstant(latt_dir));
}
}
template <Uint DIM>
Tree<DIM>::~Tree() {
for (typename std::vector<Tree<DIM> *>::iterator child= this->children.begin();
child != this->children.end(); ++child) {
delete *child;
}
}
template <Uint DIM>
void Tree<DIM>::fill(PointContainer<DIM> & points) {
// the first one always has to be split manually, because none of its corners are
// within el domaino
if (doIHaveToSplit()) {
this->split();
typename std::vector<Tree<DIM> * >::iterator child= this->children.begin();
for (; child != this->children.end() ; ++child) {
(*child)->fill(points);
}
} else {
this->create_points(points);
}
}
template <Uint DIM>
void Tree<DIM>::dump(std::string file_name) {
std::ofstream outfile;
outfile.open(file_name.c_str());
this->dump(outfile);
outfile.close();
}
template <Uint DIM>
Tree<DIM>::Tree(Real * south_west_, Uint size_, Lattice<DIM> * lattice_,
Domain<DIM> * domain_, Uint neighbours_)
:size(size_), domain(domain_), lattice(lattice_), neighbours(neighbours_) {
this->nb_boxes ++;
for (Uint i = 0; i < DIM ; ++i) {
this->south_west[i] = south_west_[i];
}
}
template <Uint DIM>
void Tree<DIM>::dump(std::ostream & stream) {
Uint nb_corners = 1<<DIM;
Real coords[DIM];
Real scaled_south_west[DIM];
Uint reorder[nb_corners];
if (DIM==2) {
reorder[0] = 0;
reorder[1] = 1;
reorder[2] = 3;
reorder[3] = 2;
} else if (DIM == 3) {
for (Uint i = 0; i < nb_corners ; ++i) {
reorder[i] = i;
}
} else {
throw(std::string("dump has been implemented for 2 and 3 D trees only"));
}
for (Uint corner_id = 0; corner_id < nb_corners ; ++corner_id) {
vector_set<DIM>(coords, this->south_west);
vector_set<DIM>(scaled_south_west, this->south_west);
for (Uint latt_dir = 0; latt_dir < DIM ; ++latt_dir) {
Uint factor = reorder[corner_id]>>latt_dir & 1;
vector_incr<DIM>(coords, this->lattice->getLatticeVector(latt_dir),
factor * this->size * this->lattice->getConstant(latt_dir));
}
for (Uint dir = 0; dir < DIM ; ++dir) {
stream << coords[dir] << "\t";
}
stream << std::endl;
}
for (Uint dir = 0; dir < DIM; ++dir) {
stream << scaled_south_west[dir] << "\t";
}
stream << std::endl << std::endl;
typename std::vector<Tree<DIM> * >::iterator child= this->children.begin();
for (; child != this->children.end() ; ++child) {
(*child)->dump(stream);
}
}
template <Uint DIM>
void Tree<DIM>::split() {
Uint new_size = size >> 1;
Uint geometric_size = new_size;
keep_max(geometric_size, static_cast<Uint>(1));
static Uint nb_children = 1<<DIM;
if (new_size == 0) {
Tree<DIM> * child = new Tree<DIM>(this->south_west, new_size, this->lattice,
this->domain, this->neighbours);
child->level = this->level+1;
this->children.push_back(child);
} else {
for (Uint child_id = 0; child_id < nb_children ; ++child_id) {
Real child_sw[DIM];
vector_set<DIM>(child_sw, this->south_west);
for (Uint latt_dir = 0; latt_dir < DIM ; ++latt_dir) {
Uint factor = child_id >> latt_dir & 1;
vector_incr<DIM>(child_sw, this->lattice->getLatticeVector(latt_dir),
factor * geometric_size*this->lattice->getConstant(latt_dir));
}
Tree<DIM> * child = new Tree<DIM>(child_sw, new_size, this->lattice,
this->domain, ~(this->neighbours | child_id));
this->children.push_back(child);
child->level = this->level+1;
}
}
}
template <Uint DIM>
bool Tree<DIM>::doIHaveToSplit() {
static Uint nb_corners = 1<<DIM;
Real coords[DIM];
if (this->size == 0 ) {
return false;
}
for (Uint corner_id = 0; corner_id < nb_corners ; ++corner_id) {
vector_set<DIM>(coords, this->south_west);
for (Uint latt_dir = 0; latt_dir < DIM ; ++latt_dir) {
Uint factor = corner_id>>latt_dir & 1;
vector_incr<DIM>(coords, this->lattice->getLatticeVector(latt_dir),
factor * this->size * this->lattice->getConstant(latt_dir));
}
Real distance_from_border = -this->domain->from_border(coords) + this->size*std::numeric_limits<Real>::epsilon();
if ( distance_from_border >=0 ) {
Real refinement = this->domain->interpolate(coords);// * this->lattice->getRepresentativeConstant();
/* got to satisfy two conditions to be split:
* 1: not too close to border
* 2: not larger than local refinement */
bool condition = this->size > refinement;
if (condition) {
return true;
}
}
}
if (this->level <= 1) {
return true;
}
return false;
}
template <Uint DIM>
void Tree<DIM>::create_points(PointContainer<DIM> & points) {
PointContainer<DIM> temp_points;
if (this->size == 0){
// in this case, the lattice needs to be inserted
this->lattice->fillAtoms(this->south_west, temp_points);
} else {
//in this case, we have to set the south_west node and maybe neighbours
temp_points.addPoint(south_west);
if (this->neighbours != 0 && false) {
Real neighbour_coords[DIM];
for (Uint neighbour = 0; neighbour < DIM ; ++neighbour) {
if (!(this->neighbours>>neighbour & 1)) {
vector_set<DIM>(neighbour_coords, this->south_west);
vector_incr<DIM>(neighbour_coords, this->lattice->getLatticeVector(neighbour),
this->size * this->lattice->getConstant(neighbour));
temp_points.addPoint(neighbour_coords);
}
}
}
}
for (Uint i = 0; i < temp_points.getSize() ; ++i) {
PointRef<DIM> current_point = temp_points.getPoint(i);
if (this->domain->inDomain(current_point)) {
points.addPoint(current_point);
}
}
}
template<Uint DIM>
void Tree<DIM>::printself(std::ostream & stream, int indent) const {
indent_space(stream, indent);
stream << "IMPLEMENT THIS";}
template<Uint DIM>
unsigned int Tree<DIM>::nb_boxes=0;
template class Tree<twoD>;
template class Tree<threeD>;
c4science · Help