tree.cpp :
Code:
#include "tetra.cpp"
template<typename T>
void print_mesh(buffer<T> & buff) {
tree<T, tetra<T>*>* printer = reinterpret_cast<tree<T, tetra<T>*>*>(buff.tree_buffer.data());
size_t distance = (buff.tree_position - buff.tree_buffer.data()) / sizeof(tree<T, tetra<T>*>);
unsigned int leaves = 0;
for (decltype(distance) i = 0; i < distance; ++i) {
if ((printer+i)->children[0] == nullptr && (printer+i)->children[1] == nullptr && (printer+i)->children[2] == nullptr && (printer+i)->children[3] == nullptr) {
++leaves;
print_tetra((printer + i)->t, "Mesh tetrahedron :");
}
}
std::cout << "\nTotal number of tetrahedra in the mesh : " << leaves << std::endl;
return;
}
template<typename T>
void print_tree(buffer<T> & buff) {
tree<T, tetra<T>*>* printer = reinterpret_cast<tree<T, tetra<T>*>*>(buff.tree_buffer.data());
size_t distance = (buff.tree_position - buff.tree_buffer.data()) / sizeof(tree<T, tetra<T>*>);
unsigned int leaves = 0;
for (decltype(distance) i = 0; i < distance; ++i) {
print_tetra((printer + i)->t, "Tree tetrahedron");
if ((printer+i)->children[0] == nullptr && (printer+i)->children[1] == nullptr && (printer+i)->children[2] == nullptr && (printer+i)->children[3] == nullptr) {
++leaves;
}
}
std::cout << "\nTotal number of leaves in the tree : " << distance << std::endl;
std::cout << "Total number of tetrahedra in the mesh : " << leaves << std::endl;
return;
}
template<typename T>
void fracture(buffer<T> & buff, std::vector<tree<T, tetra<T>*>*> & leaves, point<T>* p) {
if (debug) std::cout << "\nSize of leaf pile : " << leaves.size() << std::endl;
for (auto it = leaves.begin(); it < leaves.end(); ++it) {
tree<T, tetra<T>*>* parent = *it;
if (debug) {
print_tetra(parent->t, "Fracturing...");
}
std::array<point<T>*, 4> sample = (*it)->t->p;
point<T>* a = sample[0];
point<T>* b = sample[1];
point<T>* c = sample[2];
point<T>* d = sample[3];
std::array<std::array<point<T>*, 4> , 4> vertices = { { { a, b, c, p }, { a, c, d, p }, { a, d, b, p }, { d, c, b, p } } };
if (debug) {
std::cout << "\nThis is all the point data for a facture :\n";
for (unsigned int i = 0; i < 4; ++i) {
std::cout << "Point set index : " << i << std::endl;
print_point(vertices[i][0]);
print_point(vertices[i][1]);
print_point(vertices[i][2]);
print_point(vertices[i][3]);
std::cout << "\n";
}
}
for (unsigned int i = 0; i < 4; ++i) {
if (debug) std::cout << "Fracture index : " << i << std::endl;
if (alloc_tetra(buff, vertices[i][0], vertices[i][1], vertices[i][2], vertices[i][3])) {
continue;
}
alloc_tree(buff, reinterpret_cast<tetra<T>*>(buff.tetra_position - sizeof(tetra<T>)));
tree<T, tetra<T>*>* frac = reinterpret_cast<tree<T, tetra<T>*>*>(buff.tree_position - sizeof(tree<T, tetra<T>*>));
frac->parents[0] = parent;
if (parent)
parent->children[i] = frac;
if (debug) {
print_tetra(frac->t, "Fractured child");
std::cout << "\n";
}
}
}
return;
}
template<typename T>
void search_tree(tree<T, tetra<T>*>* root, point<T>* p, std::vector<tree<T, tetra<T>*>*> & leaves) {
std::set<tree<T, tetra<T>*>*> visited;
std::stack<tree<T, tetra<T>*>*> traversal;
traversal.push(root);
while (!traversal.empty()) {
tree<T, tetra<T>*>* curr = traversal.top();
if (!visited.insert(curr).second) {
traversal.pop();
continue;
}
if (debug) {
print_tetra( curr->t, "Visiting...");
}
if (curr->del == 0) {
if (inside_tetra(curr->t, p) == -1) {
traversal.pop();
continue;
}
unsigned int sum = 0;
for (auto i = 0; i < curr->children.size(); ++i)
if (curr->children[i])
++sum;
if (sum == 0) {
if (debug) {
print_tetra(curr->t, "Adding this to the leaf pile...");
}
leaves.push_back(curr);
traversal.pop();
continue;
} else {
//for (auto it = curr->children.end(); it < curr->children.end(); ++it)
//traversal.push(*it);
for (auto it = curr->children.crbegin(); it != curr->children.crend(); ++it)
if (*it)
traversal.push(*it);
continue;
}
} else if (curr->del == 1) {
//for (auto it = curr->children.begin(); it < curr->children.end(); ++it)
//traversal.push(*it);
for (auto it = curr->children.crbegin(); it != curr->children.crend(); ++it)
if (*it)
traversal.push(*it);
continue;
}
}
return;
}
template<typename T>
void alloc_tree(buffer<T> & buff, tetra<T>* t) {
tree<T, tetra<T>*> leaf;
for (unsigned int i = 0; i < leaf.children.size(); ++i)
leaf.children[i] = nullptr;
leaf.t = t;
memcpy(buff.tree_position, &leaf, sizeof(tree<T, tetra<T>*>));
buff.tree_position += sizeof(tree<T, tetra<T>*>);
return;
}
structures.hpp :
Code:
#include <iostream>
#include <array>
#include <vector>
#include <cstring>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <exception>
#include <string>
#include <set>
#include <stack>
/* Base structure for the graph */
template<typename T, class X>
struct tree {
std::array<struct tree<T, X>*, 4> parents;
std::array<struct tree<T, X>*, 4> children;
X t; // tetrahedron contained by node
unsigned int del;
tree(void) { del = 0; }
};
/* Point structure */
template<typename T>
struct point {
T x, y, z;
point(T a, T b, T c) : x(a), y(b), z(c) { };
};
/* Tetrahedral structure */
template<typename T>
struct tetra {
std::array<struct tetra<T>*, 4> ngb; // 4 neigbhours
std::array<struct point<T>*, 4> p; // 4 vertices
struct tree<T, struct tetra<T>*>* h; // host leaf
std::array<T, 4> a, b, c, d; // half-plane data
tetra(point<T>* a, point<T>* b, point<T>* c, point<T>* d);
};
/* Create a structure to write data into */
template<typename T>
struct buffer {
std::vector<char> point_buffer;
std::vector<char> tetra_buffer;
std::vector<char> tree_buffer;
char* point_position;
char* tetra_position;
char* tree_position;
buffer(unsigned int ppp) {
point_buffer.resize((5 + ppp) * sizeof(point<T>));
tetra_buffer.resize((4+3*10*(ppp-1))*sizeof(tetra<T>));
tree_buffer.resize((5+4*10*(ppp-1))*sizeof(tree<T, struct tetra<T>*>));
point_position = this->point_buffer.data();
tetra_position = this->tetra_buffer.data();
tree_position = this->tree_buffer.data();
}
};
/* Exception handling classes */
class co_linear : public std::exception {
public :
co_linear(std::string m = "Three points are co-linear") : msg(m) { }
~co_linear(void) throw() { }
const char* what() const throw() { return msg.c_str();}
private :
std::string msg;
};
class co_planar : public std::exception {
public :
co_planar(std::string m = "Fourth point is co-planar") : msg(m) { }
~co_planar(void) throw() { }
const char* what() const throw() { return msg.c_str();}
private :
std::string msg;
};
/* Global variables and functions */
#define U double
#define debug false
extern unsigned int num_threads;
extern unsigned int box_length;
int io(int argc, char** argv);