August 8, 2010

Today I worked some more on reformatting the fracture code. I have completed the Tetrahedron class (not tested or compiled, but the code has been transferred)

- I optimized the code in the calculating the strain/stress/derivatives (mostly eliminating redundant loops.

- Currently, the tetrahedron class is untemplated and uses floats. It also uses only House:: namespace libraries. This should be fixed in the future to take a custom data type. Also, the linear math system should either use a different library, or matrix/vector types should be passed as template parameters.

- I was thinking about making the tetrahedron class an abstract class with basic finite element functions which can be implemented with a derived class. This way, multiple element styles/numerical techniques can be tested by simply deriving with a custom element class.

- I am going to use smart pointers for the list of Node pointers which the tetrahedrons keep track of. This way, they can automatically be deleted when no tetrahedron is using them. There is an issue I need to look up before I do this.
? Can a std::vector of smart pointers resize itself when a smart pointer has no more owner? If not, then there will still be a dangling pointer in the vector which would eliminate any advantage of using smart pointers (of course, a routine per time step could go through and resize the vector).

More will be researched later.

August 7, 2010

Today I finally got my fracture stuff implemented. Currently, it is extremely slow and not optimized, but I finally implemented an approximate fracture plane. Basically, I take each node and calculate the separation tensor described in Obrien.
I think take the points within the element, and get a vector from the node to those points ( Ni - Ncurrent where Ni != Ncurrent ). I think take the dot product with the eigenvector (seperation vector) and if the absolute value is greater than 0 for 2 or more nodes in the element, then that element gets assigned to one side of the fracture plane. 5 or 6 hours later, I got all of the memory leaks worked out. Basically, after this is done, new nodes are created, old nodes may be deleted, and the node's tetrahedron index values are changed, so I have to go through the nodes and delete all the unused nodes (if they are no longer part of the mesh) and update all of the elements' index values into this vector of nodes (for each of the 4 nodes).

Also, I looked at Ogre, Openscenegraph, and Bullet. I'm going to experiment with Ogre and Openscenegraph, but Bullet seems to be a good physics library so I might see if I can extend it with the fracture stuff.

I also got all the LOP tests for work. I finally managed to fix and issue with the grid. So now it goes through each cell, takes the neighbors and uses the neighbors in its calculations.
Preliminary results: The points approximate the thin surface of each object, but they are not equally spaced out as in the paper. More experimenting will be done to see if this is an issue with my code (more than likely not) or just an artifact of using the original uniformly sampled dataset to project itself onto itself (in the paper, they use a "blanket" of regular samples to project these points causing them to not clump). The only thing I'm worried about is the clumping in general.

Tomorrow I'm going to continue to experiment with fracture. I will also read some more papers and try to flesh out sketch of thesis proposal / previous works.

Todo next: Clean up code and rewrite to be more modular

Uniform Geometry Hashed Grid

Hey all,

So recently, I was asked to design a few applications which processed large data structures such as point clouds. The algorithms which I needed to implement were mostly smoothing algorithms where for each point or data in the cloud, you had to look at all the other data around it and perform some computations with that data. Obviously, if you are using a nearest neighbor approach, this is not that bad. It still suffers some drawbacks, such as when you are not sure of the area to process on or if you simply can't use a nearest neighbors approach to reduce the computational overhead. An obvious choice would be to create some type of partitioning structure which will allow each data point to only take into account those other points within its geometric range ( an individual cell ). So when dealing with arbitrary data points where you don't know how many cells you will need or even if the data would do well with a uniform grid partition (perhaps most points lie along the same X, Y planes and don't span the Z plane that much), then a different type of grid partitioning system would need to be used in these situations. The method that I chose to use and aim to discuss here is the geometricaly hashed grid.

In this formulation, each grid cell has a coordinate set of 3 integers. When looking for a cell to add the point to, the point's position is quickly and effiently converted to the grid cell's coordinate system and added to the bucket represented by that grid cell containing all the other poitns in the cell. The grid cell's integer coordinates are used as the key in a has table linking the cell's key to a bucket containing the points within that cell. The one caveat here is that cells are added IF AND ONLY IF there are points for this cell. This means that unneeded memory (and processing) overhead are reduced by making sure that empty cells are not stored in memory. It turns out that this structure is easy to code, efficient to process, and works well with testing algorithms which require some kind of partitioning structure to run in a reasonable amount of time, but don't need a complex data structure such as an octree quite yet.

So how did I want to interface this structure? It turns out that the amount of coding for this structure can be minimized using some of C++'s STL containers/algorithms, as well as templating it for easy expansion. When I show this code, please be aware that it is not a true hash algorithm in the sense that the coordinates used to index the table or not hashed to an index value for constant time access. Instead, the hash table remains sorted in memory, so whenever a grid's cell coordinate is used to index the hash table, there is a time complexity of Log N to perform the binary search. Some hash schemes could be extended to hash the 3 integer coordinates of the grid cells into a single integer for constant time access, but I'll leave that up to you :-D

Throughout this, Point<3, CT> == a 3-dimensional point where dimensions are type CT
First is the overall structure for the grid class:

template <typename CT>
class Grid {
public:
        //
        //Constructor/Desctructors
 //

        Grid ( const std::vector< Point< 3, CT > > &, CT );
 ~Grid();

        //
 //Gets a vector a points from each cell and those within (dist)cells away from it
 //for example, if dist = 1, then all cells immediately touching the cell accessed
 //would contribute their points.  If dist = 2, then cells touching the immediate cells would contribute.
 //etc.  This is nessacary for operations needing points near the grid cells.
 //

        void GetNeighborCells( int, int, std::vector< Point< 3, CT > > & );

        // 
        //Insert single point into grid cell.
 //

        void InsertIntoGrid( const Point< 3, CT > &);

        // 
        //Retrieve a grid cell
 //
 
        std::vector< Point< 3, CT > > GetCellElements( int, int, int );

        // 
        //Gets number of cells in grid
 //

        int CurrNumberOfCells() { return mNCells; }
 
        //
        //Gets number of elements in the grid
 //

        int CurrNumberOfElements() { 
  int retVal = 0;
  for (std::map<Point<3, float>, Bucket >::const_iterator cit = mTable.begin();
   cit != mTable.end();
   ++cit)
   retVal += cit->second.data.size();
  return retVal;
 }
 
        //
        //Prints all of the index values for each grid cell (integers)
 //

        void PrintCellIndexValues(){
  for (std::map<Point<3, float>, Bucket >::const_iterator cit = mTable.begin();
   cit != mTable.end();
   ++cit)
 }
 
        //
        //Individually access each cell with a single index value.
 //

        std::vector< Point< 3, CT> > operator[](int i){
  assert ( (i+1) <= mTable.size() && (i >= 0) ); 
  //need to get iterator since map[] does not do what you would expect
  std::map< Point<3, float>, Bucket >::const_iterator cit = mTable.begin();
  for (int j = 0; j < i; j++)
   ++cit; //increment to the proper position
  std::vector < Point<3, CT> > retVal = cit->second.data;
  return (retVal);
 }
 

protected:
private:
 //Methods
 Grid ();

        //Used to hold all the point elements for each cell.
 typedef struct _Bucket{
  std::vector< Point< 3, CT > > data;
 } Bucket;

        
 //Members
 CT mEdgeLength; //Edge length of each cell
 int mNCells; //number of cells total
 std::map< Point< 3, float >, Bucket > mTable; //Actual Hash Table
};

While I'm not going to show each and ever method, I will show how the constructor sets up the initial grid, how the accessors work using the stl algorithms.

//
//Constructor
//

template <typename CT>
Grid<CT>::Grid( const std::vector< Point<3, CT> > & inPts, CT edgeLength){
 //Set # cells to 0
 mNCells = 0; //Initialize number of cells to 0
 mEdgeLength = edgeLength; //Initialize edgeLength
 Point< 3, float > tempPoint; //Used for intermediate data
 Bucket tempBucket; //Used for intermediate data
 int i, j, k; //Cell index values.

 std::pair< std::map< Point<3, float>, Bucket>::iterator, bool > check; //Temporary individual elements of the hash table.

 
 //insert the data into the cells
 if (inPts.size() == 0){
  std::cout << "Input points size == 0" << std::endl;
 }
 else{
  //Insert to see if data is in the hash table
  for (std::vector < Point< 3, CT> >::const_iterator cit = inPts.begin();
   cit != inPts.end();
   ++cit){
    
                                //Take the 3-d coordinates of the point, and 
                                //convert into cell coordinates.  This is the 
                                //simple way to go about it.
    i = floor(cit->X() / mEdgeLength);
    j = floor(cit->Y() / mEdgeLength);
    k = floor(cit->Z() / mEdgeLength);

    tempBucket.data.clear(); //clear temporary data from previous iteration.
    tempBucket.data.push_back( (*cit) ); //Add the point to the bucket
    tempPoint.X() = i;
    tempPoint.Y() = j;
    tempPoint.Z() = k;
 
                                //Since this is a map, if this insert operation fails due to
                                //The element already being in the hash table, then check.second will == false.

    check = mTable.insert( std::pair< Point<3, float>, Bucket > (tempPoint, tempBucket ) );
    //Use the return value to see if a new cell was added (previously empty)
    //or if the insertion failed and the point just needs to be added to the current
    //bucket
    if (check.second == false){ //element already in hash table,
                                //So simply add point to current bucket for cell's grid coordinate.
     (((check.first)->second).data).push_back(*cit); //Put point in bucket.
    }
    //else element successfully inserted.
    else{
     ++mNCells; //A new cell was added.
    }
  }
 }
}

And that's how the grid is set up. Basically, it converts each dimensional coordinate of each point into a 3-d grid integar coordinate. This coordinate is used to hash the grid cells. While iterating through the points, if the cells exist (the hash value already has an existing entry/bucket in the table) then the point is added to that entry's bucket. If not, the a new entry/buck are created and inserted into the std::map hash table.

To access each cell, I have created a method which will take the coordinate or index value and return a vector of points associated with that index value. Although this is simple and can be extended in many ways, it gets the job done.

//
// Get all poitns associated with hash index
//

template <typename CT>
std::vector< Point< 3, CT > > Grid<CT>::GetCellElements( int i, int j, int k ){
 Point<3, float> tempPoint;
 std::vector< Point<3, CT> > retVal;
 std::map< Point<3, float>, Bucket >::const_iterator cit;
 
 tempPoint.X() = i;
 tempPoint.Y() = j;
 tempPoint.Z() = k;
 cit = this->mTable.find( tempPoint ); //looks for value
 if ( cit != mTable.end() ){ //found the key
//  assert (cit->second.data.size() != 0); //if 0, something went horribly wrong
  retVal.insert( retVal.begin(), cit->second.data.begin(), cit->second.data.end() );
 }

 return retVal;
}

Also, the operator[] which instead takes an index value (great for iteration algorithms where the cell coordinates are not known and you simply want to iterate through cells [0] to [N] where N == number of cells in grid table) is inlined and shown up in the first code posted in this blog.

A proper destructor should be used to deallocate any memory, and if storage is a concern, this table's buckets could instead store references or pointers to an externally allocated array of points and dereference these when sorting/checking for grid coordinates. In this case, a custom sort algorithm must be provided to std::map for it to sort properly (it uses the default std::less algorithm which uses the templated classe's < operator to keep the map sorted, so make sure your point class has these operators overloaded!!!!! I hope this helps, and let me know if you have any questions or comments.

First Post

To keep with the tradition of first blog posts being short, I will in no way abandon those principles.
Hello, world!

I will post development related things on here. Tomorrow, I will post some info about the 3D hash grid using C++ and some STL concepts. Stay tuned!