No Description

engine.c 5.1KB

    /* Implementation of the Barnes-Hut algorithm */ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> #include <string.h> #include <unistd.h> #include <stdbool.h> #include "engine.h" #define THETA 0.5 void initTree(struct node_t *tree, double xmin, double xmax, double ymin, double ymax); void calcForces(struct node_t *tree, RigidMass_t *mass); //need a way to deal with masses that are on top of each other void addToTree(struct node_t *tree, RigidMass_t *masses, int n, int rlev) { //printf("n=%d\r\n", n); if(n == 2) { rlev++; } /*RigidMass_t *masses_topleft; int n_topleft = 0; RigidMass_t *masses_topright; int n_topright = 0; RigidMass_t *masses_botleft; int n_botleft = 0; RigidMass_t *masses_botright; int n_botright = 0; masses_topleft = malloc(n * sizeof(RigidMass_t)); masses_topright = malloc(n * sizeof(RigidMass_t)); masses_botleft = malloc(n * sizeof(RigidMass_t)); masses_botright = malloc(n * sizeof(RigidMass_t));*/ RigidMass_t *massQuadrants[4]; //Array of Array of massquadrants for(int i = 0; i < 4; i++) { massQuadrants[i] = malloc(sizeof(RigidMass_t) * n); } int massQuadCounts[] = {0, 0, 0, 0}; tree->n += n; double centerX = (tree->xmax + tree->xmin) / 2; double centerY = (tree->ymax + tree->ymin) / 2; // printf("Center of this node: (%f, %f)", centerX, centerY); double totalMass = 0; double weightedCenterX = 0; double weightedCenterY = 0; for(int i = 0; i < n; i++) { //add mass to tree stats totalMass += masses[i].mass; weightedCenterX += masses[i].mass * masses[i].x; weightedCenterY += masses[i].mass * masses[i].y; //categorize mass for child node double curX = masses[i].x; double curY = masses[i].y; if(curX < centerX) { //Left side if(curY < centerY) { //Top left massQuadrants[0][massQuadCounts[0]++] = masses[i]; } else { //Bot left massQuadrants[2][massQuadCounts[2]++] = masses[i]; } } else { if(curY < centerY) { //Top right massQuadrants[1][massQuadCounts[1]++] = masses[i]; } else { //Bot right massQuadrants[3][massQuadCounts[3]++] = masses[i]; } } } tree->CoM_x = weightedCenterX / totalMass; tree->CoM_y = weightedCenterY / totalMass; tree->mass = totalMass; //recursion if(rlev < 100) //can't go more than 100 levels deep. { for(int i = 0; i < 4; i++) { if(massQuadCounts[i] > 1) { tree->nodes[i] = malloc(sizeof(struct node_t)); //This gets freed when the tree is deleted double treeXmin = i % 2 == 0 ? tree->xmin : centerX; double treeXmax = i % 2 == 0 ? centerX : tree->xmax; double treeYmin = i < 2 ? tree->ymin : centerY; double treeYmax = i < 2 ? centerY : tree->ymax; initTree(tree->nodes[i], treeXmin, treeXmax, treeYmin, treeYmax); addToTree(tree->nodes[i], massQuadrants[i], massQuadCounts[i], rlev + 1); } //Free up the quadrant memory free(massQuadrants[i]); } } } void initTree(struct node_t *tree, double xmin, double xmax, double ymin, double ymax) { tree->mass = 0; tree->n = 0; tree->xmin = xmin; tree->xmax = xmax; tree->ymin = ymin; tree->ymax = ymax; tree->nodes[0] = 0; tree->nodes[1] = 0; tree->nodes[2] = 0; tree->nodes[3] = 0; } void delTree(struct node_t *tree, bool isParent) { if(tree == 0) return; //printf("delTree recursion!\r\n"); delTree(tree->nodes[0], false); delTree(tree->nodes[1], false); delTree(tree->nodes[2], false); delTree(tree->nodes[3], false); if(!isParent) { free(tree); //printf("tree freed!\r\n"); } } void barnesIterate(RigidMass_t *masses, int n, int w, int h) { printf("Attempting to make a tree...\r\n"); struct node_t tree; //quad tree encapsulates an area of 5w * 5h, centered on the screen view initTree(&tree, -2*w, 3*w, -2*h, 3*h); addToTree(&tree, masses, n, 0); printf("Tree done. Calculating forces...\r\n"); #pragma omp parallel for for(int i = 0; i < n; i++) { calcForces(&tree, &masses[i]); } printf("Deleting tree...\r\n"); delTree(&tree, true); printf("Processing tick..\r\n"); for(int i = 0; i < n; i++) { masses[i].vx += masses[i].dvx; masses[i].vy += masses[i].dvy; masses[i].dvx = 0; masses[i].dvy = 0; masses[i].x += masses[i].vx * masses[i].tick_time; masses[i].y += masses[i].vy * masses[i].tick_time; } } const double G = 0.00667408; void calcForces(struct node_t *tree, RigidMass_t *mass) { double dx = tree->CoM_x - mass->x; double dy = tree->CoM_y - mass->y; //done in many parts b/c we need dist_squared and dist when approxing anyway double dist_squared = dx*dx + dy*dy; double dist = sqrt(dist_squared); double sd = (tree->xmax - tree->xmin) / dist; bool approx = true; //if all of the children are null ptrs, we can not go any lower for(int i = 0; i < 4; i++) { if(tree->nodes[i] != 0) approx = false; } if(sd < THETA || approx) //we can approximate { //printf("s/d = %f\r\n", sd); if(dist > 1) { double v = G * tree->mass / dist_squared * mass->tick_time; mass->dvx += v * (dx) / dist; mass->dvy += v * (dy) / dist; } } else { //go down a level for(int i = 0; i < 4; i++) { if(tree->nodes[i] != 0) { calcForces(tree->nodes[i], mass); } } } }