main.cpp :
Code:
#include <iostream>
#include <vector>
#include <memory>
#include <exception>
#include "structures.h"
using std::cout;
using std::endl;
using std::vector;
using std::exception;
int tetra_count = 1;
int main(void) {
/* Allocate Cartesian distribution of particles */
vector<particle> particle_set;
for (int i=0; i<8; i++) {
for (int j=0; j<8; j++) {
for (int k=0; k<8; k++) {
particle_set.push_back(particle(i, j, k));
}
}
}
/* Create an all-encompassing root tetrahedron
and print all relevant information */
particle root_p0(0, 0, 0);
particle root_p1(21, 0, 0);
particle root_p2(0, 21, 0);
particle root_p3(0, 0, 21);
tetra* root_tetra = NULL;
try {
root_tetra = new tetra(&root_p0, &root_p1, &root_p2, &root_p3);
} catch (exception & z) {
cout << z.what() << endl;
cout << "Failed to allocate root tetrahedron" << endl;
return 1;
}
print_tetra("Root tetrahedron :", root_tetra);
int inside = 0, surface = 0, outside = 0;
for (int i=0; i<512; i++) {
int x = insideTetra(root_tetra, &particle_set[i]);
if (x == 0) inside++;
else if (x == -1) outside++;
else if (x > 0) surface++;
}
cout << "\nThere are " << inside << " points inside the root" << endl;
cout << "There are " << surface << " points on the surface of the root" << endl;
cout << "There are " << outside << " points outside the root" << endl;
cout << inside + surface << " total contained points" << endl;
/* Allocate the root of the tree */
tree *tree_root = new tree(root_tetra);
/* Begin the insertion routine, using our particle set */
vector<tree*> found;
search_tree(tree_root, &particle_set[255], found);
grow_tree(&particle_set[255], found);
cout << "Number of children the root has : " << tree_root->children.size() << endl;
/* Deallocate various references in memory */
delete(tree_root);
/* Display various information */
cout << "\nTotal tetrahedrons : " << tetra_count << endl;
cout << "\nSize of struct particle : " << sizeof(particle) << " bytes" << endl;
cout << "Size of struct tetra : " << sizeof(tetra) << " bytes" << endl;
cout << "Size of struct tree : " << sizeof(tree) << " bytes" << endl;
return 0;
}
tetra.cpp :
Code:
#include <iostream>
#include <exception>
#include <vector>
#include <string>
#include "structures.h"
using std::cout;
using std::endl;
using std::exception;
using std::vector;
using std::string;
/* Routine used to display a message and the vertices of
a tetrahedron */
void print_tetra(string message, tetra *t) {
cout << "\n" + message << endl;
for (int i=0; i<4; i++)
printf("(%.00f, %.00f, %.00f)\n", t->p[i]->x, t->p[i]->y, t->p[i]->z);
return;
}
/* Insertion routine for fracturing tetrahedra by
vertex insertion */
vector<tetra*> fracture(tetra *t, particle *p) {
vector<tetra*> fract;
if (insideTetra(t, p) == -1) {
cout << "Point not contained within tetrahedron" << endl;
return fract;
}
tetra_count--;
particle *p0 = t->p[0];
particle *p1 = t->p[1];
particle *p2 = t->p[2];
particle *p3 = t->p[3];
print_tetra("Parent tetrahedron :", t);
try {
tetra *child = new tetra(p0, p1, p2, p);
fract.push_back(child);
print_tetra("First sub-tetrahedron :", child);
tetra_count++;
} catch (exception & z) {
cout << "\n" << z.what() << endl;
cout << "Failed to construct first subtetrahedon (012)" << endl;
}
try {
tetra *child = new tetra(p0, p2, p3, p);
fract.push_back(child);
print_tetra("Second sub-tetrahedron :", child);
tetra_count++;
} catch (exception & z) {
cout << "\n" << z.what() << endl;
cout << "Failed to construct second subtetrahedon (023)" << endl;
}
try {
tetra *child = new tetra(p0, p3, p1, p);
fract.push_back(child);
print_tetra("Third sub-tetrahedron :", child);
tetra_count++;
} catch (exception & z) {
cout << "\n" << z.what() << endl;
cout << "Failed to construct third subtetrahedon (031)" << endl;
}
try {
tetra *child = new tetra(p3, p2, p1, p);
fract.push_back(child);
print_tetra("Fourth sub-tetrahedron :", child);
tetra_count++;
} catch (exception & z) {
cout << "\n" << z.what() << endl;
cout << "Failed to construct fourth subtetrahedon (321)" << endl;
}
return fract;
}
/* Routine used to determine a point's location relative to
a tetrahedron */
int insideTetra(tetra *t, particle *p) {
double d[4] = {
t->a[0]*p->x + t->b[0]*p->y + t->c[0]*p->z + t->d[0],
t->a[1]*p->x + t->b[1]*p->y + t->c[1]*p->z + t->d[1],
t->a[2]*p->x + t->b[2]*p->y + t->c[2]*p->z + t->d[2],
t->a[3]*p->x + t->b[3]*p->y + t->c[3]*p->z + t->d[3],
};
if (d[0] > 0 && d[1] > 0 && d[2] > 0 && d[3] > 0)
return 0; /* Inside */
if (d[0] < 0 || d[1] < 0 || d[2] < 0 || d[3] < 0)
return -1; /* Outside */
else {
int sum = 0;
for (int i=0; i<4; i++)
if (d[i]==0) sum++;
return sum; /* Vertex, surface or edge */
}
}
/* Constructor for each tetrahedron. We initialize the half-plane data
and adjust the orientation to be positive */
tetra::tetra(particle *p0, particle *p1, particle *p2, particle *p3) {
{ /* Plane 012 */
double x10 = p1->x - p0->x;
double y10 = p1->y - p0->y;
double z10 = p1->z - p0->z;
double x20 = p2->x - p0->x;
double y20 = p2->y - p0->y;
double z20 = p2->z - p0->z;
double a = y10*z20 - y20*z10;
double b = -(x10*z20 - x20*z10);
double c = x10*y20 - x20*y10;
if (a==0 && b==0 && c==0)
throw co_linear();
double d = a*p0->x + b*p0->y + c*p0->z;
//printf("%.00Lf*x + %.00Lf*y + %.00Lf*z = %.00Lf\n", a, b, c, d);
d = -d;
double p = a*p3->x + b*p3->y + c*p3->z + d;
if (p > 0) {
tetra::a[0] = a;
tetra::b[0] = b;
tetra::c[0] = c;
tetra::d[0] = d;
tetra::p[0] = p0;
tetra::p[1] = p1;
tetra::p[2] = p2;
tetra::p[3] = p3;
} else if (p < 0) {
tetra::a[0] = -a;
tetra::b[0] = -b;
tetra::c[0] = -c;
tetra::d[0] = -d;
tetra::p[0] = p0;
tetra::p[1] = p2;
tetra::p[2] = p1;
tetra::p[3] = p3;
} else if (p == 0)
throw co_planar();
}
{ /* Plane 023 */
double x20 = p2->x - p0->x;
double y20 = p2->y - p0->y;
double z20 = p2->z - p0->z;
double x30 = p3->x - p0->x;
double y30 = p3->y - p0->y;
double z30 = p3->z - p0->z;
double a = y20*z30 - y30*z20;
double b = -(x20*z30 - x30*z20);
double c = x20*y30 - x30*y20;
double d = a*p0->x + b*p0->y + c*p0->z;
//printf("%.00Lf*x + %.00Lf*y + %.00Lf*z = %.00Lf\n", a, b, c, d);
d = -d;
tetra::a[1] = a;
tetra::b[1] = b;
tetra::c[1] = c;
tetra::d[1] = d;
}
{ /* Plane 031 */
double x30 = p3->x - p0->x;
double y30 = p3->y - p0->y;
double z30 = p3->z - p0->z;
double x10 = p1->x - p0->x;
double y10 = p1->y - p0->y;
double z10 = p1->z - p0->z;
double a = y30*z10 - y10*z30;
double b = -(x30*z10 - x10*z30);
double c = x30*y10 - x10*y30;
double d = a*p0->x + b*p0->y + c*p0->z;
//printf("%.00Lf*x + %.00Lf*y + %.00Lf*z = %.00Lf\n", a, b, c, d);
d = -d;
tetra::a[2] = a;
tetra::b[2] = b;
tetra::c[2] = c;
tetra::d[2] = d;
}
{ /* Plane 321 */
double x23 = p2->x - p3->x;
double y23 = p2->y - p3->y;
double z23 = p2->z - p3->z;
double x13 = p1->x - p3->x;
double y13 = p1->y - p3->y;
double z13 = p1->z - p3->z;
double a = y23*z13 - y13*z23;
double b = -(x23*z13 - x13*z23);
double c = x23*y13 - x13*y23;
double d = a*p3->x + b*p3->y + c*p3->z;
//printf("%.00Lf*x + %.00Lf*y + %.00Lf*z = %.00Lf\n", a, b, c, d);
d = -d;
tetra::a[3] = a;
tetra::b[3] = b;
tetra::c[3] = c;
tetra::d[3] = d;
}
/* Don't forget to initialize the face values */
tetra::ngb[0] = NULL;
tetra::ngb[1] = NULL;
tetra::ngb[2] = NULL;
tetra::ngb[3] = NULL;
}
tree.cpp :
Code:
#include <iostream>
#include <vector>
#include "structures.h"
using std::cout;
using std::endl;
using std::vector;
/* Basic search algorithm for our tree */
void search_tree(tree *root, particle *p, vector<tree*> & stack) {
if (insideTetra(root->t, p) == -1)
return;
if (root->children.empty() && root->occ == 0) {
root->occ = 1;
stack.push_back(root);
return;
}
for (vector<tree*>::iterator it = root->children.begin(); it < root->children.end(); it++) {
if (insideTetra((*it)->t, p) != -1)
search_tree(*it, p, stack);
}
return;
}
/* Fracture all collected leaves and link them by
shared faces */
void grow_tree(particle *p, vector<tree*> & stack) {
vector<tetra*> total_fract;
for (vector<tree*>::iterator it = stack.begin(); it < stack.end(); it++) {
vector<tetra*> fract = fracture((*it)->t, p);
for (vector<tetra*>::iterator it2 = fract.begin(); it2 < fract.end(); it2++) {
(*it)->children.push_back(new tree(*it2));
total_fract.push_back(*it2);
}
for (int i=0; i<4; i++)
if ((*it)->t->ngb[i])
total_fract.push_back((*it)->t->ngb[i]);
}
cout << "Number of fractures in total_fract " << total_fract.size() << endl;
inside_sort(total_fract);
return;
}
neighbour.cpp :
Code:
#include <iostream>
#include <vector>
#include <cmath>
#include "structures.h"
using std::cout;
using std::endl;
using std::vector;
void big_switch(tetra *t1, tetra *t2, int hash) {
switch (hash) {
/* 012 block */
case 27 : //012 - 012
t1->ngb[0] = t2;
t2->ngb[0] = t1;
break;
case 243 : //012 - 023
t1->ngb[0] = t2;
t2->ngb[1] = t1;
break;
case 81 : //012 - 031
t1->ngb[0] = t2;
t2->ngb[2] = t1;
break;
case 729 : //012 - 321
t1->ngb[0] = t2;
t2->ngb[3] = t1;
break;
/* 023 block */
case 125 : //023 - 012
t1->ngb[1] = t2;
t2->ngb[0] = t1;
break;
case 3125 : //023 - 023
t1->ngb[1] = t2;
t2->ngb[1] = t1;
break;
case 625 : //023 - 031
t1->ngb[1] = t2;
t2->ngb[2] = t1;
break;
case 15625 : //023 - 321
t1->ngb[1] = t2;
t2->ngb[3] = t1;
break;
/* 031 block */
case 64 : //031 - 012
t1->ngb[2] = t2;
t2->ngb[0] = t1;
break;
case 1024 : //031 - 023
t1->ngb[2] = t2;
t2->ngb[1] = t1;
break;
case 256 : //031 - 031
t1->ngb[2] = t2;
t2->ngb[2] = t1;
break;
case 4096 : //031 - 321
t1->ngb[2] = t2;
t2->ngb[3] = t1;
break;
/* 321 block */
case 216 : //321 - 012
t1->ngb[3] = t2;
t2->ngb[0] = t1;
break;
case 7776 : //321 - 023
t1->ngb[3] = t2;
t2->ngb[1] = t1;
break;
case 1296 : //321 - 031
t1->ngb[3] = t2;
t2->ngb[2] = t1;
break;
case 46656 : //321 - 321
t1->ngb[3] = t2;
t2->ngb[3] = t1;
break;
}
return;
}
void inside_sort(vector<tetra*> & fract) {
for (vector<tetra*>::iterator it = fract.begin(); it < fract.end()-1; it++) {
for (vector<tetra*>::iterator it2 = it+1; it2<fract.end(); it2++) {
int sum1 = 0, sum2 = 0, count = 0;
for (int i=0; i<4; i++) {
for (int j=0; j<4; j++) {
if ((*it)->p[i] == (*it2)->p[j]) {
sum1 += i;
sum2 += j;
count++;
}
}
}
if (count == 3) {
big_switch(*it, *it2, pow(sum1, sum2));
}
}
}
for (vector<tetra*>::iterator it = fract.begin(); it < fract.end(); it++) {
print_tetra("The following are the neighbours of ", *it);
for (int i=0; i<4; i++) {
if ((*it)->ngb[i]) {
cout << "Neighbouring index : " << i << endl;
print_tetra("Neighbouring tetrahedron", (*it)->ngb[i]);
}
}
}
return;
}
structures.h :
Code:
#include <iostream>
#include <string>
#include <exception>
#include <vector>
using std::cout;
using std::endl;
using std::string;
using std::exception;
using std::vector;
/* Basic particle structure */
struct particle {
public :
double x;
double y;
double z;
particle(double a, double b, double c) : x(a), y(b), z(c) { }
};
/* Basic tetrahedron structure */
struct tetra {
public :
struct tetra *ngb[4];
struct particle *p[4];
double a[4], b[4], c[4], d[4];
tetra(particle *p0, particle *p1, particle *p2, particle *p3);
};
/* Node structure for our tree */
struct tree {
public :
vector<struct tree*> children;
struct tetra *t;
unsigned occ : 1;
tree(tetra *insert) {
t = insert;
occ = 0;
}
~tree(void) {
delete(t);
}
};
/* Exception handling classes */
class co_linear : public exception {
public :
co_linear(string m = "Three points are co-linear") : msg(m) { }
~co_linear(void) throw() { }
const char* what() const throw() { return msg.c_str();}
private :
string msg;
};
class co_planar : public exception {
public :
co_planar(string m = "Fourth point is co-planar") : msg(m) { }
~co_planar(void) throw() { }
const char* what() const throw() { return msg.c_str();}
private :
string msg;
};
/* Various global variables and functions */
extern int tetra_count;
int insideTetra(tetra *t, particle *p);
vector<tetra*> fracture(tetra *t, particle *p);
void print_tetra(string message, tetra *t);
void search_tree(tree *root, particle *p, vector<tree*> & stack);
void grow_tree(particle *p, vector<tree*> & stack);
void inside_sort(vector<tetra*> & fract);
Makefile :
Code:
MKFILE = Makefile
GCC = g++ -g -O3 -Wall -std=c++11 -Wextra -lpthread
CSOURCE = main.cpp tetra.cpp tree.cpp neighbour.cpp
CHEADER = structures.h
OBJECTS = ${CSOURCE:.cpp=.o}
EXECBIN = regulus
all : ${EXECBIN}
${EXECBIN} : ${OBJECTS}
${GCC} -o $@ ${OBJECTS} -lm
%.o : %.cpp
${GCC} -c $<
clean :
- rm ${OBJECTS} ${EXECBIN}
again :
${MAKE} clean all