You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

224 lines
5.1 KiB

  1. /* Implementation of the Barnes-Hut algorithm */
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include <time.h>
  6. #include <string.h>
  7. #include <unistd.h>
  8. #include <stdbool.h>
  9. #include "engine.h"
  10. #define THETA 0.5
  11. void initTree(struct node_t *tree, double xmin, double xmax, double ymin, double ymax);
  12. void calcForces(struct node_t *tree, RigidMass_t *mass);
  13. //need a way to deal with masses that are on top of each other
  14. void addToTree(struct node_t *tree, RigidMass_t *masses, int n, int rlev)
  15. {
  16. //printf("n=%d\r\n", n);
  17. if(n == 2)
  18. {
  19. rlev++;
  20. }
  21. /*RigidMass_t *masses_topleft;
  22. int n_topleft = 0;
  23. RigidMass_t *masses_topright;
  24. int n_topright = 0;
  25. RigidMass_t *masses_botleft;
  26. int n_botleft = 0;
  27. RigidMass_t *masses_botright;
  28. int n_botright = 0;
  29. masses_topleft = malloc(n * sizeof(RigidMass_t));
  30. masses_topright = malloc(n * sizeof(RigidMass_t));
  31. masses_botleft = malloc(n * sizeof(RigidMass_t));
  32. masses_botright = malloc(n * sizeof(RigidMass_t));*/
  33. RigidMass_t *massQuadrants[4]; //Array of Array of massquadrants
  34. for(int i = 0; i < 4; i++)
  35. {
  36. massQuadrants[i] = malloc(sizeof(RigidMass_t) * n);
  37. }
  38. int massQuadCounts[] = {0, 0, 0, 0};
  39. tree->n += n;
  40. double centerX = (tree->xmax + tree->xmin) / 2;
  41. double centerY = (tree->ymax + tree->ymin) / 2;
  42. // printf("Center of this node: (%f, %f)", centerX, centerY);
  43. double totalMass = 0;
  44. double weightedCenterX = 0;
  45. double weightedCenterY = 0;
  46. for(int i = 0; i < n; i++)
  47. {
  48. //add mass to tree stats
  49. totalMass += masses[i].mass;
  50. weightedCenterX += masses[i].mass * masses[i].x;
  51. weightedCenterY += masses[i].mass * masses[i].y;
  52. //categorize mass for child node
  53. double curX = masses[i].x;
  54. double curY = masses[i].y;
  55. if(curX < centerX)
  56. {
  57. //Left side
  58. if(curY < centerY)
  59. {
  60. //Top left
  61. massQuadrants[0][massQuadCounts[0]++] = masses[i];
  62. }
  63. else
  64. {
  65. //Bot left
  66. massQuadrants[2][massQuadCounts[2]++] = masses[i];
  67. }
  68. }
  69. else
  70. {
  71. if(curY < centerY)
  72. {
  73. //Top right
  74. massQuadrants[1][massQuadCounts[1]++] = masses[i];
  75. }
  76. else
  77. {
  78. //Bot right
  79. massQuadrants[3][massQuadCounts[3]++] = masses[i];
  80. }
  81. }
  82. }
  83. tree->CoM_x = weightedCenterX / totalMass;
  84. tree->CoM_y = weightedCenterY / totalMass;
  85. tree->mass = totalMass;
  86. //recursion
  87. if(rlev < 100) //can't go more than 100 levels deep.
  88. {
  89. for(int i = 0; i < 4; i++)
  90. {
  91. if(massQuadCounts[i] > 1)
  92. {
  93. tree->nodes[i] = malloc(sizeof(struct node_t)); //This gets freed when the tree is deleted
  94. double treeXmin = i % 2 == 0 ? tree->xmin : centerX;
  95. double treeXmax = i % 2 == 0 ? centerX : tree->xmax;
  96. double treeYmin = i < 2 ? tree->ymin : centerY;
  97. double treeYmax = i < 2 ? centerY : tree->ymax;
  98. initTree(tree->nodes[i], treeXmin, treeXmax, treeYmin, treeYmax);
  99. addToTree(tree->nodes[i], massQuadrants[i], massQuadCounts[i], rlev + 1);
  100. }
  101. //Free up the quadrant memory
  102. free(massQuadrants[i]);
  103. }
  104. }
  105. }
  106. void initTree(struct node_t *tree, double xmin, double xmax, double ymin, double ymax)
  107. {
  108. tree->mass = 0;
  109. tree->n = 0;
  110. tree->xmin = xmin;
  111. tree->xmax = xmax;
  112. tree->ymin = ymin;
  113. tree->ymax = ymax;
  114. tree->nodes[0] = 0;
  115. tree->nodes[1] = 0;
  116. tree->nodes[2] = 0;
  117. tree->nodes[3] = 0;
  118. }
  119. void delTree(struct node_t *tree, bool isParent)
  120. {
  121. if(tree == 0) return;
  122. //printf("delTree recursion!\r\n");
  123. delTree(tree->nodes[0], false);
  124. delTree(tree->nodes[1], false);
  125. delTree(tree->nodes[2], false);
  126. delTree(tree->nodes[3], false);
  127. if(!isParent)
  128. {
  129. free(tree);
  130. //printf("tree freed!\r\n");
  131. }
  132. }
  133. void barnesIterate(RigidMass_t *masses, int n, int w, int h)
  134. {
  135. printf("Attempting to make a tree...\r\n");
  136. struct node_t tree;
  137. //quad tree encapsulates an area of 5w * 5h, centered on the screen view
  138. initTree(&tree, -2*w, 3*w, -2*h, 3*h);
  139. addToTree(&tree, masses, n, 0);
  140. printf("Tree done. Calculating forces...\r\n");
  141. #pragma omp parallel for
  142. for(int i = 0; i < n; i++)
  143. {
  144. calcForces(&tree, &masses[i]);
  145. }
  146. printf("Deleting tree...\r\n");
  147. delTree(&tree, true);
  148. printf("Processing tick..\r\n");
  149. for(int i = 0; i < n; i++)
  150. {
  151. masses[i].vx += masses[i].dvx;
  152. masses[i].vy += masses[i].dvy;
  153. masses[i].dvx = 0;
  154. masses[i].dvy = 0;
  155. masses[i].x += masses[i].vx * masses[i].tick_time;
  156. masses[i].y += masses[i].vy * masses[i].tick_time;
  157. }
  158. }
  159. const double G = 0.00667408;
  160. void calcForces(struct node_t *tree, RigidMass_t *mass)
  161. {
  162. double dx = tree->CoM_x - mass->x;
  163. double dy = tree->CoM_y - mass->y;
  164. //done in many parts b/c we need dist_squared and dist when approxing anyway
  165. double dist_squared = dx*dx + dy*dy;
  166. double dist = sqrt(dist_squared);
  167. double sd = (tree->xmax - tree->xmin) / dist;
  168. bool approx = true;
  169. //if all of the children are null ptrs, we can not go any lower
  170. for(int i = 0; i < 4; i++)
  171. {
  172. if(tree->nodes[i] != 0) approx = false;
  173. }
  174. if(sd < THETA || approx) //we can approximate
  175. {
  176. //printf("s/d = %f\r\n", sd);
  177. if(dist > 1)
  178. {
  179. double v = G * tree->mass / dist_squared * mass->tick_time;
  180. mass->dvx += v * (dx) / dist;
  181. mass->dvy += v * (dy) / dist;
  182. }
  183. }
  184. else
  185. {
  186. //go down a level
  187. for(int i = 0; i < 4; i++)
  188. {
  189. if(tree->nodes[i] != 0)
  190. {
  191. calcForces(tree->nodes[i], mass);
  192. }
  193. }
  194. }
  195. }