3 * @purpose Implementation of Barnes Hut algorithms
4 * @author Sam Moore (20503628) 2012
13 unsigned nested_used = 0;
18 * @function Node_SetSpace
19 * @purpose Set the position and space occupied by a node
20 * @param node - The Node
21 * @param x - The vertex with the smallest coordinate values
22 * @param size - The dimensions of the cube that the Node represents
24 void Node_SetSpace(Node * node, float x[DIMENSIONS], float size[DIMENSIONS])
26 for (unsigned i = 0; i < DIMENSIONS; ++i)
29 node->size[i] = size[i];
34 * @function Node_Destroy
35 * @purpose Confuse Java programmers
36 * *ahem* Cleanup the memory used by the Node and its children
37 * @param node - The Node
39 void Node_Destroy(Node * node)
45 //printf("Destroy node %p\n", (void*)node);
47 if (node->body != NULL)
56 Prepare_Threads(SUB_DIVISIONS);
57 #pragma omp parallel for
59 for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
61 Node_Destroy(node->children[i]);
62 free(node->children[i]);
63 node->children[i] = NULL;
69 * @function Node_FindBodies
70 * @purpose Search for bodies in a system within a Node, and add them to the Node's array of bodies.
73 * @returns The total number of bodies found from s within n
75 unsigned Node_FindBodies(Node * n, System * s)
81 // reset the body array first
83 //Prepare_Threads(options.nested_threads);
84 //#pragma omp parallel for
86 for (unsigned i = 0; i < n->allocated; ++i)
89 n->N = 0; // Don't forget this!
91 if (n->parent == NULL)
94 //Prepare_Threads(options.nested_threads);
95 //#pragma omp parallel for
97 for (unsigned a = 0; a < s->N; ++a)
99 if (Node_ContainsBody(n, s->body+a))
100 Node_AddBody(n, s->body+a);
106 //Prepare_Threads(options.nested_threads);
107 //#pragma omp parallel for
109 for (unsigned a = 0; a < n->parent->allocated; ++a)
111 if (Node_ContainsBody(n, n->parent->body[a]))
112 Node_AddBody(n, n->parent->body[a]);
115 //printf("Node %p contains %d bodies\n", (void*)(n), n->N);
120 * @function Node_Subdivide
121 * @purpose Subdivide a Node into more nodes and fill them with bodies
122 * @param node - The Node
123 * @param s - The system which the bodies are from
124 * @returns The number of bodies contained in all children of node after this is finished
126 unsigned Node_Subdivide(Node * node, System * s)
128 if (Node_FindBodies(node, s) <= 1) // First find all bodies in the node
129 return node->N; // If there was one or less, no need to subdivide
131 unsigned count = 0; // stores the return value of the function
132 float size[DIMENSIONS]; // size of the children of this node
133 for (unsigned i = 0; i < DIMENSIONS; ++i)
134 size[i] = node->size[i] / 2.0;
137 // Prepare_Threads(SUB_DIVISIONS);
138 // #pragma omp parallel for
140 for (unsigned a = 0; a < SUB_DIVISIONS; ++a)
142 if (node->children[a] == NULL)
143 node->children[a] = Node_Create(node); //Create child node if necessary
147 // Determine the placement of each child node in space
148 // Annoyingly long switch statement.
157 x[0] = node->x[0] + size[0];
163 x[1] = node->x[1] + size[1];
169 x[2] = node->x[2] + size[2];
172 x[0] = node->x[0] + size[0];
173 x[1] = node->x[1] + size[1];
177 x[0] = node->x[0] + size[0];
179 x[2] = node->x[2] + size[2];
184 x[1] = node->x[1] + size[1];
185 x[2] = node->x[2] + size[2];
188 x[0] = node->x[0] + size[0];
189 x[1] = node->x[1] + size[1];
190 x[2] = node->x[2] + size[2];
194 Node_SetSpace(node->children[a], x, size); // Setup the child in space
195 unsigned sub = Node_Subdivide(node->children[a], s); // Recurse for the child
203 * @function Node_AddBody
204 * @purpose Terrify third year computer science students
205 * *ahem* Add a Body* to the node's dynamic Body* array. Dynamically. Got that?
206 * I love pointers! They are fun!
207 * @param node - The Node; what did you think it was?
208 * @param b - The Body to add
209 * @returns - The index in the array that the Body* now occupies
211 unsigned Node_AddBody(Node * node, Body * b)
214 for (i = 0; i < node->allocated; ++i) // Search for an empty spot in the array
216 if (node->body[i] == NULL) // Got one!
218 node->body[i] = b; // Fill it!
224 // No empty spots found :(
225 // The index got outside the array's allocated space, so allocate some more!
226 node->allocated = 2 * (node->allocated + 1); // Doubling the space each time, like what C++ std::vector does
228 if (node->body == NULL)
229 node->body = (Body**)(calloc(node->allocated, sizeof(Body*))); // First Body* added, create new array
230 else // Subsequent allocations use realloc
232 // realloc is expensive. That's why I double the space each time, instead of just adding a constant.
233 // Of course, in the worst case, we are wasting half the memory after the last realloc
234 // But performance wise it is better.
235 node->body = (Body**)(realloc(node->body, sizeof(Body*) * node->allocated)); // Resize existing array
236 memset(node->body+i, 0, sizeof(Body*) * (node->allocated - i)); // set all the new elements to NULL
239 //We now have a valid position in the array that is not filled
240 node->body[i] = b; // So fill it
241 node->N += 1; // don't forget to update this
242 return i; //index in array of Body*
247 * @function Generate_Tree
248 * @purpose Generate the tree from a System
249 * @param s - The System
250 * @param n - The root node
251 * @returns The number of bodies in the node (should be the same as in the system)
253 unsigned Generate_Tree(System * s, Node * n)
255 for (unsigned i = 0; i < DIMENSIONS; ++i) // set initial size of node
261 Node_FitSystem(n, s); // size node so that it contains all bodies in the system
264 Node_Subdivide(n, s);
270 * @function Node_FitSystem
271 * @purpose Resize a Node until all bodies in a System are within it
272 * @param n - The Node
273 * @param s - The System
275 void Node_FitSystem(Node * n, System * s)
277 // I'm sure there is a faster way of doing this
278 // But I like this algorithm.
281 while (Node_FindBodies(n, s) < s->N) // Find bodies in the node...
283 // As long as the number found is less than the number in the system...
285 for (unsigned i = 0; i < DIMENSIONS; ++i)
288 n->size[i] *= 2.0; // Double the size of the Node!
294 * @function Update_Tree
295 * @purpose Update the tree (it must already have been Generated with Generate_Tree)
296 * @param s - The System
297 * @param n - The root node of the tree
298 * @returns The number of bodies in the tree
300 unsigned Update_Tree(System * s, Node * n)
302 Node_FitSystem(n, s);
304 Node_Subdivide(n, s);
309 * @function Node_CalculateMass
310 * @purpose Calculate the mass and centre of mass of a node
311 * @param node - The node
312 * @returns The mass of the node
314 float Node_CalculateMass(Node * node)
317 if (node->N <= 0) // No bodies in the node
320 if (node->N == 1) // Only one body in the node
322 Body * b = Node_Body(node, 0);
323 node->mass = b->mass;
324 for (unsigned i = 0; i < DIMENSIONS; ++i)
325 node->com[i] = b->x[i];
326 //printf("Body %p has mass %f, center of mass at %f %f %f\n", (void*)(node), node->mass, node->com[0], node->com[1], node->com[2]);
331 // For more than 1 body in the node, calculate mass and com from child nodes
333 for (unsigned i = 0; i < DIMENSIONS; ++i) // init com
338 Prepare_Threads(SUB_DIVISIONS);
339 #pragma omp parallel for
341 for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
343 if (node->children[i] != NULL)
345 float m = Node_CalculateMass(node->children[i]); // Recurse for children
347 for (unsigned j = 0; j < DIMENSIONS; ++j)
348 node->com[j] += m * node->children[i]->com[j];
351 for (unsigned i = 0; i < DIMENSIONS; ++i)
352 node->com[i] = node->com[i] / node->mass; // weighted average for centre of masses
354 //printf("Node %p has mass %f, center of mass at %f %f %f\n", (void*)(node), node->mass, node->com[0], node->com[1], node->com[2]);
359 * @function Node_Forces
360 * @purpose Calculate the forces on a body due to a Node
361 * @param a - The Body
362 * @param n - The node
363 * @param theta - The approximation factor; 0.0 for no approximation.
364 * @returns The depth reached by the recursive algorithm
366 unsigned Node_Forces(Body * a, Node * n, float theta)
372 if (n->N == 1) // external node with one body in it
374 Body * b = Node_Body(n, 0);
375 if (b == a) //ignore self
379 Body_Force(a, b); // calculate force
383 // internal node with more than one body in it...
389 float s = n->size[0] + n->size[1] + n->size[2] / 3.0;
391 for (unsigned i = 0; i < DIMENSIONS; ++i)
392 d += square(a->x[i] - n->com[i]);
395 if (s / d <= theta) // approximation is "valid"
397 // Calculate force based on centre of mass of the node
398 con = G * a->mass * n->mass / square(d);
400 for (unsigned i = 0; i < DIMENSIONS; ++i)
401 a->F[i] += gd * (n->com[i] - a->x[i]);
405 // approximation "invalid" because node is too close
406 // Then recursively calculate forces due to the node's children
410 Prepare_Threads(SUB_DIVISIONS);
411 #pragma omp parallel for
413 for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
415 if (n->children[i] != NULL)
417 count += Node_Forces(a, n->children[i], theta);
424 * @function Node_Create
425 * @purpose Create a Node
426 * @param parent - The parent of the Node (NULL if the node is the root of the tree)
427 * @returns The new Node
429 Node * Node_Create(Node * parent)
431 Node * result = (Node*)(calloc(1, sizeof(Node)));
433 result->parent = parent;
435 /* // In theory, using calloc instead of malloc should make everything zero
436 // So I don't need to do it manually
437 // Kept code here anyway. Because I don't trust gcc.
440 result->allocated = 0;
442 for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
443 result->children[i] = NULL;
451 * @function Node_Body
452 * @purpose Identify a Body within a Node by its position in the array of Body*
453 * Because elements of the array may be NULL in general, n->body[index] may fail
454 * Actually I think in my code it will work in all cases, because I never set elements of the array to NULL.
455 * It's not worth the risk anyway, I might change my code.
457 * @param n - The Node
458 * @param index - The position of the Body* in the array
459 * @returns the Body* from the array
461 Body * Node_Body(Node * n, unsigned index)
464 Body * result = NULL;
466 // Can't parallelise this, due to the break statement
467 for (unsigned i = 0; i < n->allocated; ++i) // Go through the array...
469 if (n->body[i] != NULL)
470 count += 1; // Count each non-empty entry
473 result = n->body[i]; // return element once counter expires
482 * @function Draw_Node
483 * @purpose Draw a node as a series of rectangular prisms
484 * Draws all children of the node too. But only nodes with N == 1 will be shown.
485 * @param n - The Node
487 void Draw_Node(Node * n)
493 glColor3f(0.0f, 1.0f, 0.0f);
494 glVertex3f(n->x[0] + n->size[0], n->x[1] + n->size[1], n->x[2]);
495 glVertex3f(n->x[0], n->x[1] + n->size[1], n->x[2]);
496 glVertex3f(n->x[0], n->x[1] + n->size[1], n->x[2] + n->size[2]);
497 glVertex3f(n->x[0]+n->size[0], n->x[1]+n->size[1], n->x[2]+n->size[2]);
499 glVertex3f(n->x[0] + n->size[0], n->x[1], n->x[2] + n->size[2]);
500 glVertex3f(n->x[0], n->x[1] , n->x[2] + n->size[2]);
501 glVertex3f(n->x[0], n->x[1], n->x[2]);
502 glVertex3f(n->x[0] + n->size[0], n->x[1], n->x[2]);
504 glVertex3f(n->x[0]+n->size[0], n->x[1]+n->size[1], n->x[2]+n->size[2]);
505 glVertex3f(n->x[0], n->x[1]+n->size[1], n->x[2]+n->size[2]);
506 glVertex3f(n->x[0], n->x[1], n->x[2]+n->size[2]);
507 glVertex3f(n->x[0]+n->size[0], n->x[1], n->x[2]+n->size[2]);
509 glVertex3f(n->x[0]+n->size[0], n->x[1], n->x[2]);
510 glVertex3f(n->x[0], n->x[1], n->x[2]);
511 glVertex3f(n->x[0], n->x[1]+n->size[1], n->x[2]);
512 glVertex3f(n->x[0]+n->size[0], n->x[1]+n->size[1], n->x[2]);
514 glVertex3f(n->x[0]+n->size[0], n->x[1]+n->size[1], n->x[2]);
515 glVertex3f(n->x[0]+n->size[0], n->x[1]+n->size[1], n->x[2] + n->size[2]);
516 glVertex3f(n->x[0]+n->size[0], n->x[1], n->x[2] + n->size[2]);
517 glVertex3f(n->x[0]+n->size[0], n->x[1], n->x[2]);
520 //printf("Node %p contains %u bodies in volume %f\n", (void*)n, n->N, n->size[0] * n->size[1] * n->size[2]);
521 for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
523 if (n->children[i] != NULL && n->children[i]->N > 0)
524 Draw_Node(n->children[i]);
530 * @function Node_ContainsBody
531 * @purpose Indicate whether body is within a node
532 * @param n - The Node
533 * @param b - The Body
534 * @returns true iff the body was in the node
536 bool Node_ContainsBody(Node * n, Body * b)
540 for (unsigned i = 0; i < DIMENSIONS; ++i)
542 float p = b->x[i] - n->x[i];
543 if (p < 0.00 || p > n->size[i])
552 * @function Prepare_Threads
553 * @purpose Tell OMP to spawn up to the specified number of threads on the next parallel section
554 * @param max - Spawn at most this many threads; in reality, less will be spawned
556 void Prepare_Threads(unsigned max)
560 if (options.nested_threads == 1)
562 omp_set_num_threads(1);
565 //printf("Thread %d in nest %d thinking of spawning more threads...\n", omp_get_thread_num(), omp_get_level());
566 //printf("There are %d nested threads, %d in current level\n", nested_used, omp_get_num_threads());
567 //printf("Got past check\n");
569 if (max > options.nested_threads)
570 max = options.nested_threads;
572 max = options.nested_threads;
577 unsigned available = options.nested_threads - nested_used;
583 omp_set_num_threads(max);
585 omp_set_lock(&nested_lock);
586 nested_used += max-1;
587 omp_unset_lock(&nested_lock);