webber@porthos.rutgers.edu (Bob Webber) (04/04/88)
In article <609@imagine.PAWL.RPI.EDU>, kyriazis@pawl20.pawl.rpi.edu (George Kyriazis) writes: > > Hello world. I am about to start a project about simulating rigid > bodies or fluids made out of particles, and I want to ask some questions > before beginning. ... Also calculating closest > neighborhood in 2-d or even 3-d for the next time step from the previous > neighborhood is difficult. Are there any data structures or algorithms > that can (easily) handle this kind of thing? The simplest approach is to store your points in a three-dimensional array, i.e., in essence bucket sort them. This leads to the question of how big should the buckets be. A general rule of them is that a bucket shouldn't contain more than 20 items (a sloppy generalization of the rule of them for when it is worthwhile to keep items sorted). Counts on the number of particles in a particle system mentioned in Particle Systems - A Technique for Modeling a Class of Fuzzy Objects [William T. Reeves, ACM Transactions on Computer Graphics, April 83, 91-108] ranged from 25,000 particles to 750,000 particles. In Approximate and Probabilistic Algorithms for Shading and Rendering Structured Particle Systems [William T. Reeves, ACM SIGGRAPH'85, 313-322.] the counts go up to 1.1 million particles. Depending on the size of the particle system you choose to work on and how fast the particles move per simulation step and the effective size of the main memory on the machine you do your work on and the size of the disk on the machine you do your work on, the bucket sort will either: 1) work fine or 2) you will find yourself thrashing alot or 3) you will find that this organization just wastes too much space. 2 and 3 would both be helped if you find an organization that wastes less space. 2 would also be helped if you found an organization that had better page fault performance. To address the last point first, assuming that you process the array in the same order as your compiler lays it out and that you ignore the effects of cells more than a few cell widths away and that your particles don't move more than a cell width or so in a given time step, you probably can't do much better with cell layout. However, if you find yourself obliged to make more aggressive probes into the structure, it is worthwhile remembering that it has been repeatedly observed that the linear ordering output by the function Traverse(Array,Width,0,0) defined as: Traverse(Array,Width,x,y) <= if (Width == 1) output Array[x][y] else begin Traverse(Array,Width/2,x,y) Traverse(Array,Width/2,x+(Width/2),y) Traverse(Array,Width/2,x,y+(Width/2)) Traverse(Array,Width/2,x+(Width/2),y+(Width/2)) end places cells that are geometrically near each other closer to each other in the linear addrress space of computer memory than the linear ordering output by the function RowColumn defined as: RowColumn(Array,Width) <= for i ranging from 0 to (Width - 1) for j ranging from 0 to (Width - 1) output Array[x][y] This result extends to three-dimensional arrays as well. Thus by laying your array out in ``Traverse'' order, you can expect fewer page faults due to probes into cells of moderate distance away. If the three-dimensional array is too verbose a data structure to use, then the next simplest is a hierarchically organized three-dimensional array usually referred to as an octree. An introduction to the various concepts surrounding the usage of such a data structure can be found in The Quadtree and Related Hierarchical Data Structures [Hanan Samet, Computing Surveys, June 1984, 187-260.] as well as in the upcoming May and July issues of the IEEE Computer Graphics & Applications journal. Briefly, in two-dimensions, you could imagine the cell Array[X][Y] (where X and Y are n-bit unsigned integers denoted X1.X2.X3...Xn and Y1.Y2.Y3...Yn) as being assigned the key Y1.X1.Y2.X2.Y3.X3...Yn.Xn that is 2n-bits long. The standard unsigned numeric ordering of such keys corresponds to the preferred geometric ordering of array cells mentioned above. Now, let us adopt the policy that if the total number of particles in the cells with the keys Y1.X1.Y2.X2...Ym.Xm.0.0, Y1.X1.Y2.X2...Ym.Xm.0.1, Y1.X1.Y2.X2...Ym.Xm.1.0, and Y1.X1.Y2.X2...Ym.Xm.1.1 is less than 10 then we will place them in an entry with the 2m-bit key Y1.X1.Y2.X2...Ym.Xm. [Naturally the four 2*(m+1) bit keys and entries are removed from the collection.] Apply this policy recursively until no more reductions can be performed. Storing all of these key/entry pairs in a B-tree [any data structures textbook, e.g., The Art of Computer Programming - Volume 3: Sorting and Searching by Donald E. Knuth, 1973] sorted by the numeric order of the key works works quite well. The obvious extensions of this idea to three dimensions would be what you want (sometimes referred to as a linear octree). Finally, there is the question of actually implementing nearest-neighbor. One can implement a true nearest neighbor search or one can implement a heuristic approximation of nearest neighbor (see Techniques for Reducing Pen Plotting Time by David P. Anderson, ACM Transactions on Computer Graphics, July 1983, 197-212) or one can simply find all points within distance D of your particle and work with them. In order to catch more global particle interactions, I assume you will be using a particle hierarchy as discussed in the Reeves articles cited above. Depending on the nature of your simulation, you may or may not be able to use that hierarchy to further partition the problem (cutting down on the number of points that you need to process at any given interaction calculation of any given time step). It should all be quite interesting. Keep us up to date on your progress. ---- BOB (webber@athos.rutgers.edu ; rutgers!athos.rutgers.edu!webber)