Thread: Let's save a little face... (tetrahedral mesh algorithm help)

Hello All,

I'm trying to allocate a tetrahedral mesh through the use of a tree. We have a root tetrahedron and vertices are inserted into it which fractures it. In fact, every insertion fractures each tetrahedron it is inserted in. So for a direct insertion, 1 tetrahedron becomes 4 because a new tetrahedron is formed out of the inserted index and each respective face of the parent.

Right now, let's just worry about direct insertions and we call this a 1-to-4 flip. But there is a such thing as a n-to-2n flip (multiple edges are shared) so my code was written to work in a more generic fashion.

I'm going to post all the code I have so you all have complete context and hopefully it won't be too much lol. But it's long so if you're scared of long code, I'd leave this thread now with a return value of 1.

What my code does, it searches the tree for valid leafs. Each leaf is a valid tetrahedron in the mesh. We use the tree for point location (which tetrahedron contains the vertex we wish to insert) so I gather up all the valid leaves first and then I fracture each gathered leaf.

I then take these fractured children and connect their faces and it's this part of the algorithm where I'm wondering if you guys can help me overcome the O(n^2) limitations of it. As of now, I'm just using vectors and iterators and it seems to be working well for what I have currently but I want to see if I can make it faster.

Here's useful information to help you understand my code.

We have a tetrahedron. So we have 4 vertices and 4 faces. We denote the faces in this order, 0 - 012, 1 - 023, 2 - 031, 3 - 321 where face 0 is made out of vertices 0, 1 and 2. Face 1 is made out of vertices 0, 2 and 3 and so on and so forth.

Right now, I'm using a direct comparison method between vertices of two differing tetrahedra but the orders aren't the same across tetrahedra even if they share the same vertices so I use a summation method. Notice that 0+1+2 = 3, 0+2+3 = 5, 0+3+1 = 4, 3+2+1 = 6 and none of the values equal each other. So the crux of my algorithm is to take the sum of the respective indices of each shared vertex and use that to connect faces.

So like tetra_1->ngb[0] = tetra_2 and tetra_2->ngb[3] = tetra_1 is how this would look. All I really wanna do is find out which faces two tetrahedra share.

I'm going to post the code in the next post to make it easier to read but please, if something isn't clear then please ask questions because I'd love to explain it as I'd love to try and get some algorithmic improvements here or just improvements for my code in general.

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

And the link map for a 1-to-4 flip like what we have here was done by me from hand a couple of months ago and it looks like this,

We denote the new tetrahedra with 1, 2, 3 and 4. And 1.012 denotes the 012 face of tetrahedron 1.

Code:

1.012 = NULL;
1.023 = 2.031;
1.031 = 3.023;
1.321 = 4.321;

2.012 = NULL;
2.023 = 3.031;
2.031 = 1.023;
2.321 = 4.031;

3.012 = NULL;
3.023 = 1.031;
3.031 = 2.023;
3.321 = 4.023;

4.012 = NULL;
4.023 = 3.321;
4.031 = 2.321;
4.321 = 1.321;

So yeah, that's how the full 1-to-4 mapping works. Notice that the 012 face of each child is NULL because that's an externally facing face so if the parent has neighbours, they will be inherited but at the moment we're just worrying about the root itself being split into 4.