course/semester2/pprog/assignment1/nbody-bh/tree.c~

   1 /**
   2  * @file "tree.c"
   3  * @purpose Implementation of Barnes Hut algorithms
   4  * @author Sam Moore (20503628) 2012
   5  */
   6
   7 #include "tree.h"
   8 #include "math.h"
   9 #include <string.h>
  10 #include "graphics.h"
  11
  12 #ifdef USE_THREADS
  13 unsigned nested_used = 0;
  14 #endif //USE_THREADS
  15
  16
  17 /**
  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
  23  */
  24 void Node_SetSpace(Node * node, float x[DIMENSIONS], float size[DIMENSIONS])
  25 {
  26         for (unsigned i = 0; i < DIMENSIONS; ++i)
  27         {
  28                 node->x[i] = x[i];
  29                 node->size[i] = size[i];
  30         }
  31 }
  32
  33 /**
  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
  38  */
  39 void Node_Destroy(Node * node)
  40 {
  41
  42         if (node == NULL)
  43                 return;
  44
  45         //printf("Destroy node %p\n", (void*)node);
  46
  47         if (node->body != NULL)
  48         {
  49                 free(node->body);
  50                 node->body = NULL;
  51         }
  52         node->allocated = 0;
  53         node->N = 0;
  54
  55         #ifdef USE_THREADS
  56                 Prepare_Threads(SUB_DIVISIONS);
  57         #pragma omp parallel for
  58         #endif //USE_THREADS
  59         for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
  60         {
  61                 Node_Destroy(node->children[i]);
  62                 free(node->children[i]);
  63                 node->children[i] = NULL;
  64         }
  65
  66 }
  67
  68 /**
  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.
  71  * @param n The Node
  72  * @param s The System
  73  * @returns The total number of bodies found from s within n
  74  */
  75 unsigned Node_FindBodies(Node * n, System * s)
  76 {
  77
  78         if (n->body != NULL)
  79         {
  80
  81                 // reset the body array first
  82                 #ifdef USE_THREADS
  83                 //Prepare_Threads(options.nested_threads);
  84                 //#pragma omp parallel for
  85                 #endif //USE_THREADS
  86                 for (unsigned i = 0; i < n->allocated; ++i)
  87                         n->body[i] = NULL;
  88         }
  89         n->N = 0; // Don't forget this!
  90
  91         if (n->parent == NULL)
  92         {
  93                 #ifdef USE_THREADS
  94                 //Prepare_Threads(options.nested_threads);
  95                 //#pragma omp parallel for
  96                 #endif //USE_THREADS
  97                 for (unsigned a = 0; a < s->N; ++a)
  98                 {
  99                         if (Node_ContainsBody(n, s->body+a))
 100                                 Node_AddBody(n, s->body+a);
 101                 }
 102         }
 103         else
 104         {
 105                 #ifdef USE_THREADS
 106                 //Prepare_Threads(options.nested_threads);
 107                 //#pragma omp parallel for
 108                 #endif //USE_THREADS
 109                 for (unsigned a = 0; a < n->parent->allocated; ++a)
 110                 {
 111                         if (Node_ContainsBody(n, n->parent->body[a]))
 112                                 Node_AddBody(n, n->parent->body[a]);
 113                 }
 114         }
 115         //printf("Node %p contains %d bodies\n", (void*)(n), n->N);
 116         return n->N;
 117 }
 118
 119 /**
 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
 125  */
 126 unsigned Node_Subdivide(Node * node, System * s)
 127 {
 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
 130
 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;
 135
 136         #ifdef USE_THREADS
 137         //      Prepare_Threads(SUB_DIVISIONS);
 138         // #pragma omp parallel for
 139         #endif //USE_THREADS
 140         for (unsigned a = 0; a < SUB_DIVISIONS; ++a)
 141         {
 142                 if (node->children[a] == NULL)
 143                         node->children[a] = Node_Create(node);  //Create child node if necessary
 144
 145                 float x[DIMENSIONS];
 146
 147                 // Determine the placement of each child node in space
 148                 // Annoyingly long switch statement.
 149                 switch (a)
 150                 {
 151                         case 0:
 152                                 x[0] = node->x[0];
 153                                 x[1] = node->x[1];
 154                                 x[2] = node->x[2];
 155                                 break;
 156                         case 1:
 157                                 x[0] = node->x[0] + size[0];
 158                                 x[1] = node->x[1];
 159                                 x[2] = node->x[2];
 160                                 break;
 161                         case 2:
 162                                 x[0] = node->x[0];
 163                                 x[1] = node->x[1] + size[1];
 164                                 x[2] = node->x[2];
 165                                 break;
 166                         case 3:
 167                                 x[0] = node->x[0];
 168                                 x[1] = node->x[1];
 169                                 x[2] = node->x[2] + size[2];
 170                                 break;
 171                         case 4:
 172                                 x[0] = node->x[0] + size[0];
 173                                 x[1] = node->x[1] + size[1];
 174                                 x[2] = node->x[2];
 175                                 break;
 176                         case 5:
 177                                 x[0] = node->x[0] + size[0];
 178                                 x[1] = node->x[1];
 179                                 x[2] = node->x[2] + size[2];
 180
 181                                 break;
 182                         case 6:
 183                                 x[0] = node->x[0];
 184                                 x[1] = node->x[1] + size[1];
 185                                 x[2] = node->x[2] + size[2];
 186                                 break;
 187                         case 7:
 188                                 x[0] = node->x[0] + size[0];
 189                                 x[1] = node->x[1] + size[1];
 190                                 x[2] = node->x[2] + size[2];
 191                                 break;
 192                 }
 193
 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
 196                 count += sub;
 197         }
 198         return count;
 199 }
 200
 201
 202 /**
 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
 210  */
 211 unsigned Node_AddBody(Node * node, Body * b)
 212 {
 213         unsigned i = 0;
 214         for (i = 0; i < node->allocated; ++i) // Search for an empty spot in the array
 215         {
 216                 if (node->body[i] == NULL) // Got one!
 217                 {
 218                         node->body[i] = b; // Fill it!
 219                         node->N += 1;
 220                         return i; // return!
 221                 }
 222         }
 223
 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
 227
 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
 231         {
 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
 237         }
 238
 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*
 243
 244 }
 245
 246 /**
 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)
 252  */
 253 unsigned Generate_Tree(System * s, Node * n)
 254 {
 255         for (unsigned i = 0; i < DIMENSIONS; ++i) // set initial size of node
 256         {
 257                 n->x[i] = -1.0;
 258                 n->size[i] = 2.0;
 259         }
 260
 261         Node_FitSystem(n, s); // size node so that it contains all bodies in the system
 262
 263
 264         Node_Subdivide(n, s);
 265
 266         return n->N;
 267 }
 268
 269 /**
 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
 274  */
 275 void Node_FitSystem(Node * n, System * s)
 276 {
 277         // I'm sure there is a faster way of doing this
 278         // But I like this algorithm.
 279
 280
 281         while (Node_FindBodies(n, s) < s->N) // Find bodies in the node...
 282         {
 283                 // As long as the number found is less than the number in the system...
 284
 285                 for (unsigned i = 0; i < DIMENSIONS; ++i)
 286                 {
 287                         n->x[i] *= 2.0;
 288                         n->size[i] *= 2.0; // Double the size of the Node!
 289                 }
 290         }
 291 }
 292
 293 /**
 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
 299  */
 300 unsigned Update_Tree(System * s, Node * n)
 301 {
 302         Node_FitSystem(n, s);
 303
 304         Node_Subdivide(n, s);
 305         return n->N;
 306 }
 307
 308 /**
 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
 313  */
 314 float Node_CalculateMass(Node * node)
 315 {
 316         node->mass = 0.0;
 317         if (node->N <= 0) // No bodies in the node
 318                 return 0.0;
 319
 320         if (node->N == 1) // Only one body in the node
 321         {
 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]);
 327                 return node->mass;
 328         }
 329
 330
 331         // For more than 1 body in the node, calculate mass and com from child nodes
 332
 333         for (unsigned i = 0; i < DIMENSIONS; ++i) // init com
 334                 node->com[i] = 0.0;
 335
 336
 337         #ifdef USE_THREADS
 338                 Prepare_Threads(SUB_DIVISIONS);
 339         #pragma omp parallel for
 340         #endif //USE_THREADS
 341         for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
 342         {
 343                 if (node->children[i] != NULL)
 344                 {
 345                         float m = Node_CalculateMass(node->children[i]); // Recurse for children
 346                         node->mass += m;
 347                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 348                                 node->com[j] += m * node->children[i]->com[j];
 349                 }
 350         }
 351         for (unsigned i = 0; i < DIMENSIONS; ++i)
 352                 node->com[i] = node->com[i] / node->mass; // weighted average for centre of masses
 353
 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]);
 355         return node->mass;
 356 }
 357
 358 /**
 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
 365  */
 366 unsigned Node_Forces(Body * a, Node * n, float theta)
 367 {
 368         if (n->N <= 0)
 369                 return 0;
 370
 371
 372         if (n->N == 1) // external node with one body in it
 373         {
 374                 Body * b = Node_Body(n, 0);
 375                 if (b == a) //ignore self
 376                 {
 377                         return 0;
 378                 }
 379                 Body_Force(a, b); // calculate force
 380                 return 1;
 381         }
 382
 383         // internal node with more than one body in it...
 384
 385
 386         float con;
 387         float gd;
 388         float d = 0;
 389         float s = n->size[0] + n->size[1] + n->size[2] / 3.0;
 390
 391         for (unsigned i = 0; i < DIMENSIONS; ++i)
 392                 d += square(a->x[i] - n->com[i]);
 393         d = sqrt(d);
 394
 395         if (s / d <= theta) // approximation is "valid"
 396         {
 397                 // Calculate force based on centre of mass of the node
 398                 con = G * a->mass * n->mass / square(d);
 399                 gd = con / d;
 400                 for (unsigned i = 0; i < DIMENSIONS; ++i)
 401                         a->F[i] += gd * (n->com[i] - a->x[i]);
 402                 return 1;
 403         }
 404
 405         // approximation "invalid" because node is too close
 406         // Then recursively calculate forces due to the node's children
 407
 408         unsigned count = 0;
 409         #ifdef USE_THREADS
 410                 Prepare_Threads(SUB_DIVISIONS);
 411         #pragma omp parallel for
 412         #endif //USE_THREADS
 413         for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
 414         {
 415                 if (n->children[i] != NULL)
 416                 {
 417                         count += Node_Forces(a, n->children[i], theta);
 418                 }
 419         }
 420         return count;
 421 }
 422
 423 /**
 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
 428  */
 429 Node * Node_Create(Node * parent)
 430 {
 431         Node * result = (Node*)(calloc(1, sizeof(Node)));
 432
 433         result->parent = parent;
 434
 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.
 438
 439         result->N = 0;
 440         result->allocated = 0;
 441         result->body = 0;
 442         for (unsigned i = 0; i < SUB_DIVISIONS; ++i)
 443                 result->children[i] = NULL;
 444         */
 445
 446         return result;
 447
 448 }
 449
 450 /**
 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.
 456  *
 457  * @param n - The Node
 458  * @param index - The position of the Body* in the array
 459  * @returns the Body* from the array
 460  */
 461 Body * Node_Body(Node * n, unsigned index)
 462 {
 463         unsigned count = 0;
 464         Body * result = NULL;
 465
 466         // Can't parallelise this, due to the break statement
 467         for (unsigned i = 0; i < n->allocated; ++i) // Go through the array...
 468         {
 469                 if (n->body[i] != NULL)
 470                         count += 1; // Count each non-empty entry
 471                 if (count > index)
 472                 {
 473                         result = n->body[i]; // return element once counter expires
 474                         break;
 475                 }
 476         }
 477
 478         return result;
 479 }
 480
 481 /**
 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
 486  */
 487 void Draw_Node(Node * n)
 488 {
 489
 490         if (n->N == 1)
 491         {
 492         glBegin(GL_QUADS);
 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]);
 498
 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]);
 503
 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]);
 508
 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]);
 513
 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]);
 518         glEnd();
 519         }
 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)
 522         {
 523                 if (n->children[i] != NULL && n->children[i]->N > 0)
 524                         Draw_Node(n->children[i]);
 525         }
 526
 527 }
 528
 529 /**
 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
 535  */
 536 bool Node_ContainsBody(Node * n, Body * b)
 537 {
 538         if (b == NULL)
 539                 return false;
 540         for (unsigned i = 0; i < DIMENSIONS; ++i)
 541         {
 542                 float p = b->x[i] - n->x[i];
 543                 if (p < 0.00 || p > n->size[i])
 544                         return false;
 545         }
 546         return true;
 547 }
 548
 549 #ifdef USE_THREADS
 550 /**
 551  * @function Prepare_Threads
 552  * @purpose Tell OMP to spawn up to the specified number of threads on the next parallel section
 553  * @param max - Spawn at most this many threads; in reality, less will be spawned
 554  */
 555 void Prepare_Threads(unsigned max)
 556 {
 557         //return;
 558
 559         if (options.nested_threads == 1)
 560         {
 561                 omp_set_num_threads(1);
 562                 return;
 563         }
 564         //printf("Thread %d in nest %d thinking of spawning more threads...\n", omp_get_thread_num(), omp_get_level());
 565         //printf("There are %d nested threads, %d in current level\n", nested_used, omp_get_num_threads());
 566         //printf("Got past check\n");
 567
 568         if (max > options.nested_threads)
 569                 max = options.nested_threads;
 570         if (max <= 0)
 571                 max = options.nested_threads;
 572
 573
 574
 575
 576         unsigned available = options.nested_threads - nested_used;
 577         if (available <= 0)
 578                 available = 1;
 579         if (available < max)
 580                 max = available;
 581
 582         omp_set_num_threads(max);
 583
 584         omp_set_lock(&nested_lock);
 585         nested_used += max-1;
 586         omp_unset_lock(&nested_lock);
 587
 588 }
 589 #endif //USE_THREADS